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TITLE OF THE INVENTION 

SYSTEM AND METHOD OF WAVEFRONT SENSING 

CROSS-REFERENCE TO RELATED APPLICATIONS 

This application claims priority under 35 U.S.C. § 1 19 to co-pending U.S. Provisional 
Patent Application No. 60/394,232, filed on 9 July 2002 in the name of Daniel M. Topa, the 
entirety of which is hereby incorporated by reference herein for all purposes as if fully set forth 
herein. 

BACKGROUND AND SUMMARY OF THE INVENTION 
[0001] Technical Field. 

[0002] This invention pertains to the field of wavefront sensing methods and devices, and more 
particularly, to a method and system for more precisely determining the location of a focal spot in 
a wavefront sensor, for more precisely determining the location of vertices or boundaries of 
lenslets in a lenslet array with respect to pixels of a wavefront sensor, for using this information 
to reconstruct an arbitrary wavefront with the wavefront sensor, and for more efficiently 
performing the calculations to perform the reconstruction. 
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[0003] Description. 

[0004] A light wavefront may be defined as the virtual surface delimited by all possible rays 
having an equal optical path length from a spatially coherent source. For example, the wavefront 
of light emanating from a point light source is a sphere (or a partial sphere where light from the 
point light source is limited to emission along a small range of angles). Meanwhile, the 
wavefront created by a collimating lens mounted at a location one focal length away from a point 
source is a plane. A wavefront may be planar, spherical, or have some arbitrary shape dictated by 
other elements of an optical system through which the light is transmitted or reflected. 
[0005] A number of different wavefront sensors and associated methods are known. Among 
these are the Shack-Hartmann wavefront sensor, which will now be described in greater detail. 
[0006] FIG. 1 shows a basic configuration of a Shack-Hartmann wavefront sensor 180. The 
Shack-Hartmann wavefront sensor 180 comprises a micro-optic lenslet array 182 that breaks an 
incoming beam into multiple focal spots 188 falling on an optical detector 184. Typically, the 
optical detector 184 comprises a pixel array, for example, a charge-coupled device (CCD) 
camera. The lenslets of the lenslet array 182 dissect an incoming wavefront and creates a pattern 
of focal spots on a charge-couple device (CCD) array. The lenslet array typically contains 
hundreds or thousands of lenslets, each on the size scale of a hundred microns. Most Shack- 
Hartmann sensors are assembled such that the CCD array is in the focal plane of the lenslet array. 
[0007] A Shack-Hartmann wavefront sensor uses the fact that light travels in a straight line, to 
measure the wavefront of light. By sensing the positions of the focal spots 188, the propagation 
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vector of the sampled light can be calculated for each lenslet of the lenslet array 182. The 
wavefront can be reconstructed from these vectors. 

[0008] To better understand one or more aspects of this invention, it is worthwhile to discuss 
the operation of a Shack-Hartmann wavefront sensor in more detail. 

[0009] In typical operation, a reference beam (e.g., a plane wave) is first imaged onto the array 
and the location of the focal spots is recorded. Then, the wavefront of interest is imaged onto the 
array, and a second set of focal spot locations is recorded. FIGs. 2A-F illustrate this process. 
[00010] Usually, some optical system delivers a wavefront onto the lenslet array which samples 
the wavefront over the tiny regions of each lenslet. The lenslets should be much smaller than the 
wavefront variation, that is, the wave-front should be isoplanatic over the sampled region. When 
the CCD array is in the focal plane of the lenslet array, each lenslet will create a focal spot on the 
CCD array. The location of these focal spots is the critical measurement, for this reveals the 
average of the wavefront slopes across each region. That is, the shift in the focal spot is 
proportional to the average of the slopes of the wavefront over the region sampled by the lenslet. 
Software may compute the shift in each focal spot. 

[00011] If the wavefront is not isoplanatic, the quality of the focal spot erodes rapidly and it 
becomes more difficult to determine the peak location. 

[00012] However, where the isoplanatic condition is satisfied and where the focal spot shift is 
consistent with the small angle approximation of Fresnel, then the focal spot shift is exactly 
proportional to the average of the wavefront slope over the lenslet. 
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[00013] The incident wavefront is then reconstructed from the measurements of the average of 
the slopes for the hundreds or thousands of lenslets in the array. 

[00014] It was stated above that the shift in the focal spot is proportional to the average of the 
slopes of the wavefront over the region sampled by the lenslet. This will now be shown as 
follows: 

[00015] We begin by examining the types of wavefronts incident upon a lenslet and calculating 
their spatial irradiance distribution in the focal plane. We then resolve what this distribution tells 
us about the average value of the slopes of the sampled wavefront. 

[00016] As noted above, most Shack-Hartmann sensors are assembled such that the CCD array 
is in the focal plane of the lenslet array. Therefore the appropriate complex amplitude 
distribution, Ff(u, v), of the field in the focal plane is given by the Fraunhofer diffraction pattern: 

1) F f (u,v)= — j j^(x,j/)exp[^ — {xu+ yv)j p(x,y)dxdy 

where the pupil function describes the lens aperture. The function is 1 inside the lens aperture 
and 0 outside of the aperture. In the example that follows we will build the pupil function using 
the rectangle function, which is defined as: 

0 for \x\ ) \ 
2) H*)=\\ M \x\ = j 

1 for \x\ < \ 
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[00017] The field F^u, v)incident upon the lenslet is related to the wavefront v|/ (x 5 y) incident 
upon the lenslet through: 

3) Ft(x,y)=exp(-W(kx 9 ky)) 

[00018] where k is the wavenumber, 2n/X. The spatial irradiance distribution in the focal plane 
is related to the field amplitude by: 

[00019] The first step was to model the portion of the wavefront sampled by the lenslet. We 
begin with planar wavefronts and advance to wavefronts with more structure. Consider the first 
six simple cases in a Taylor series: 

Vi 0 { x >y)= x 

V 20 (x,y)=x 2 

[00020] These six wavefronts were then analyzed using equation (2). Before presenting the 

results for the corresponding intensity distributions, a few functions must be defined. First is the 

sampling function: 

, 0 sin* 
(5) smcx = 
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[00021] then the error function: 



(6) 



2 x 

erf{x)^- r \e t2 dt 



[00022] and the imaginary error function: 

(7) erfi(x) = - ierf (a) 
[00023] and finally the exponential integral function: 

OO -I 2 

(8) ei(x)=~l^—dt 

-2 t 

[00024] The intensity distributions can now be written in terms of these general functions. The 
first is the well-known result: 



(9) 



i\A x >y) 



fx) 



sine 



smc 



( nd_ ^ 



[00025] The next two results are also for planar wavefronts with tilt. They are: 

J nd \ 



(10) /,„(*,>)» 



sine 



y 



and 



(11) I n (x,y) = 



,2Y 



sine 



sine 



f nd / 



(y-f) 



[00026] We expect the symmetry under interchange of x and y in equations (10) and (11) since 
space is isotropic and the lenslet has a discrete symmetry under rotations by integer multiples of 
90°. 
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[00027] For the case of nonplanar wavefronts the results quickly become more complicated. To 
present these results, first define the intermediate variables: 

2zdf+Xy n 2ndf-Xy 

(12) a = — ' and 8= — - 

2fX H 2fX 

[00028] The intensity distribution can now be written as: 

(13) I 2o{x>y)=-l[j) {erfi(za)+ erfi(zj3))(erfi(iza) + erfi(izfi))sm C 2 [^-x } 

[00029] For both l2o(x,y) and Io2(x,y) we see the symmetry in the solutions under interchange of 
x and y. The intermediate variables now become: 

2ndf + Xx , 2ndf - Xx 

04) a -—ipr md ^—yr~ 

and the irradiance distribution in the focal plane is: 

(15) Ia(x,y)= (erfi(za)+ erft(zj3))(erfi(iza) + erfi(iz/]))smc 2 ^ y 

[00030] Finally we consider the case where the input wavefront exhibits some torsion. First 
define the intermediate variables: 

litdf + Xx , Ijtdf - Xx 

(16) a ^~Tfx— ™ d ai = ~2]x~ 

a 2zdf+Xy 2ndf - Xy 

(17) &—yx— md ^ = ~W~ 

[00031] Then the irradiance distribution can be written as: 
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(18) I u (x,y)=[~^j (cf^J-^^Aj-^-^J+^-w^)) 

x (ei(-ia t fi i2 )- ei(ia^ 2 )- d(ia 2 fi i )+ ei(-iaj 2 )) 
[00032] Clearly the solutions in equations (9)-(l 1) are even functions and they are separable. 
That is, the irradiance distributions for these three functions can be written X(x)Y(y). The next 
three functions are not separable as the function arguments mix the variables x and y. However, a 
careful analysis of the constituent functions in equations (6)-(8) reveals that the combinations in 
equations (13), (15) and (18) are even functions also. 

[00033] Now we address the issue of how the amplitude of these input wavefronts affects the 
focal spot position. Consider that the wavefront has an amplitude k. The incident wavefronts 
become: 

V u (x,y)=Kxy 

[00034] These new wavefronts were again propagated through the lenslet using equation (1). 
As one might expect, there was no change in Ioo(x,y). The irradiance distributions for the odd 
parity terms become: 
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(19) 



I M>y) = [ : r 



sinc 



^(*-*/))sinc 2 



nd 



y 



and 



(20) 



sine 



[00035] For the higher order cases, there is a k slight rewrite. For example for I 2 o(x), the 
intermediate variables become: 



(21) 



a = 



IndKf + Xy 



and 



IndKf - Ay 



2fA """ r 2fA 

and the irradiance distribution in equation (13) is then multiplied by k" 1/2 . For the torsion case y 
i i(x, y) the intermediate variables are: 



(22) 



(23) 



ar,= 



fix' 



IndKf + Xx 
~lfX 

IndKf + Xy 



and 



and 



a 2 = 



IndKf - /tat 
IndKf - Ay 



2fX ' " 2/1 

and the irradiance distribution in equation 18 is then multiplied by k" 2 . 

[00036] However, the only peaks that shift position are Iio(x,y) and Ioi(x,y). And both of these 
shifts are exactly in the focal plane. The question now becomes how does this compare to the 
average value of the wavefront slopes? 
[00037] The average value of the wavefront slope is given by: 



d/2 d/2 



J \d x \if(x,y)dxdy 



(24) 



m x = 



-d/2 -d/2 
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and similarly for m y . The six input wavefronts and their average slopes m = (m x , m y ) are: 

¥\o{*>y)= kk m T = (0,0) 

Voi( x >y)= *y m T = (K,o) 

V2o( x >y)= ™ 2 m T =(0 9 fc) 
V n (x,y)= Kxy m r = (o,0) 

Vm(x>y)= *y 2 m T =(0,6) 
[00038] So for these six cases, which are the most physically relevant, the average of the 
wavefront slope exactly describes the shift of the peak in the focal plane accounting for the 
distance offset of the focal length. Because of this exact relationship, it is convenient for us to 
define the focal spot location as being the location of the peak intensity of the focal spot. 
[00039] For the purposes of this discussion, we define "isoplanatic" as the condition where the 
wavefront is well approximated over an area the size of the lenslet, by a plane wave. 
[00040] This analysis buttresses the lore of Shack-Hartmann wavefront sensing. The shift in the 
focal spot position is directly proportional to the average of the wavefront slope across the 
lenslet. In the small angle limit inherent in the Fresnel approximation, this is an exact result for 
the isoplanatic cases and also true for the higher terms examined, although focal spot location is 
ambiguous in these cases as well as for additional higher order terms. An important issue that has 
not been discussed is how combinations of these terms behave. This is an important issue that is 
quite relevant physically. 
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[00041] FIGs. 3A-F and 4A-F show some exemplary incident wavefronts and output irradiance 
distributions in the focal plane. The lenslet is a square of size d = 280 \x and has a focal length f = 
28 mm; the wavelength X = 325 nm. 

[00042] FIGS. 3A-C show the functional form of some exemplary incident wavefronts which all 
have a generally planar structure. FIGs. 3D-F graphically shows the corresponding wavefronts 
incident upon the lenslet. Note that in these examples over 95% of the light is confined to an 
area the size of a lenslet. The x and y axes are in microns. The vertical scale in the first column is 
also in microns. The vertical axis on the second column represents intensity and is in arbitrary 
units. 

[00043] FIGs. 4A-C show a series of increasingly more complicated incident wavefronts. 
Clearly the focal spots are degrading rapidly where the incident wavefront has a non-planar 
structure. Notice the concept of a peak is ambiguous here and a center of mass computation of 
focal spot location could be highly problematic. Also, only a small portion of the light falls 
within the lenslet area. The x and y axes are in microns. The vertical scale in the first column is 
also in microns. The vertical axis on the second column represents intensity and is in arbitrary 
units. 

[00044] FIGs. 3A-F and 4A-F validate the well known rule of thumb: the wavefront sensor 
works best in the isoplanatic limit. In other words, the lenslet size should be much smaller than 
the variations one is attempting to measure. Of course practical wavefront sensor design includes 
many other considerations too and one cannot shrink the lenslet size arbitrarily. But the message 
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is that in the isoplanatic limit we expect to find well-defined focal spots that the center of mass 
can handle with higher precision. 

[00045] So, to recap, a Shack-Hartrnann wavefront sensor is a powerful tool for optical analysis, 
however the precision of the wavefront data that it produces may be less than is desired at times, 
and the speed and ease by which an incident wavefront may be reconstructed are also sometimes 
less than desirable. 

[00046] First, it is difficult to precisely locate the focal spots produced by the lenslet array. 
Although the center of mass may be computed as an approximation of the focal spot location, 
more precision is desired. The basic problem with center of mass algorithms lies with deciding 
which pixels to include in the calculation. All cameras have a background level of electronic 
noise. The center of mass algorithm is especially sensitive to noise away from the brightest pixel 
because of weighting of the calculation by distance. To exclude this effect, a threshold is usually 
applied that excludes all pixels below some absolute brightness level. This threshold may be 
further refined after this preliminary step according to how bright the light is under a lenslet. To 
see the effect of this threshold, a calculation can be done using the center of mass. The sinc- 
squared intensity pattern can be overlaid with a discrete section of CCD pixels. The continuous 
intensity pattern is broken into discrete pixelated values as counts. Then the center of mass 
algorithm can be applied to the pixelated values. Then a calculation can be done where the sinc- 
squared pattern is translated into small steps across the fixed CCD. This corresponds to the 
situation where tilt is added in small steps to the beam as it hits the lenslet. The center of mass 
location can be plotted against the known tilt added to the beam. When this process is done, it 
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will be seen that the center of mass will follow a straight line up to the point where a new pixel 
starts to receive enough light to be counted (or another pixel got too little light.). At that point, 
there will be an abrupt jump in the line. The appearance will be a line with lots of little kinks in 
it. On the average, the input translation versus center of mass calculation will produce a line 
with a slope of one. So it is evident that if there are kinks in the line, the slopes of all the little 
straight sections of lines must be something different than the ideal slope of one. The kinks and 
the non-unity piecewise slopes indicates that center of mass calculations have a fundamental flaw 
that makes them unsuitable for being the basis of wavefront measurements if extremely high 
accuracies are desired. 

[00047] Secondly, in general, a lenslet boundary will not line up neatly with a pixel boundary. 
Instead the boundary line will likely cross through pixels. For precise wavefront measurements, 
it is important to precisely determine the location of the vertices or boundaries of lenslets in a 
lenslet array with respect to pixels (e.g., CCD array) of a corresponding wavefront sensor. This 
knowledge is needed in order to assign locations of the boundaries of the lenslets to locations on 
the CCD. 

[00048] Third, it is difficult to take a series of hundreds or thousands of measurements of the 
average of the slopes and reconstitute the incident wavefront. This process of taking the 
measurements and reconstituting the wavefront incident upon the lenslet array is called 
reconstruction and can be quite computationally intensive and time-consuming. 
[00049] Accordingly, it would be desirable to provide a system and method of more accurately 
locating the focal spots in a wavefront sensor. It would also be desirable to provide a system and 
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method for to precisely determine the location of the vertices or boundaries of lenslets in a lenslet 
array with respect to pixels (e.g., CCD array) of a corresponding wavefront sensor. Furthermore, 
it would be desirable to provide a system and method for reconstructing a wavefront from 
measurement data provided by a pixel array in a wavefront sensor. Even furthermore, it would 
be desirable to provide a system and method for efficiently computing polynomials for defining a 
wavefront detected by a wavefront sensor. And, it would be desirable to provide such a system 
and method which overcomes one or more disadvantages of the prior art. 
[00050] The present invention comprises a system and method for sensing a wavefront. 
[00051] In one aspect of the invention, a method of determining a location of a focal spot on a 
detector array improves upon a "pure" center of mass calculation. When a near plane wave 
illuminates a square lens, the focal spot has a sinc-squared shape. This knowledge can be 
exploited to generate a number of improvements to calculating the focal spot peak location. 
[00052] Several techniques show improvement over pure center of mass methods for locating a 
focal spot. One technique uses the knowledge that the sinc-squared function has minimum in it 
before the side lobes appear. The focal spot location algorithm can be written so these minima 
are located, and then the pixels that are identified as being proper to include in the center of mass 
calculation will be those that lie within the minima. Another technique is that the shape of the 
sinc-square function can be fit to values that show up in each pixel. This can be done with least 
squares method or by "sliding" methods. Then the equation for the sinc-squared gives the 
location of the peak. Still another refinement is to use a different weighting than a linear weight 
from center in the calculation of center of mass. The effect is to reduce weight given to dimmer 
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pixels. In particular, this weight can be matched to provide piecewise linear slopes of one, 
although this reduces the dynamic range of the sensor unless another method is found to 
eliminate the kinks. 

[00053] In another aspect of the invention, a method of determining a location and size of a 
lenslet with respect to a detector array exploits a-priori knowledge of the spacing of the lenslets 
in a lenslet array. Such a-priori knowledge can be produced by direct measurement of the lenslet 
arrays by microphotography, by knowledge of the photolithography process, or by other means. 
In a particularly advantageous manufacturing process, the spacing between lenslet is extremely 
regular. Accordingly, when one compares algorithms that calculate focal spot locations, the 
better algorithm is the one that comes up with a the locations in the most regular grid. This can 
be done on a theoretical basis using the sinc-squared patterns as inputs, or on experimental data 
using a plane wave to illuminate the lenslet array. Further refinements can be obtained by adding 
tilt to the plane wave and doing least squares fit to determine spacings. 

[00054] In a similar way, non-uniform grids can be defined if it is desired to manufacture lenslet 
arrays with non-uniform spacings. 

[00055] In yet another aspect of the invention, a method of reconstructing a wavefront from a 
plurality of focal spots produced on a detector array, includes. 

[00056] In still another aspect of the invention, a method of determining a set of polynomials 
defining a wavefront from a plurality of focal spots produced on a detector array, includes. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[00057] FIG. 1 shows a basic configuration of a Shack-Hartmann wavefront sensor. 
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[00058] FIGs. 2A-F illustrate a reference beam and a wavefront of interest being imaged onto a 
detector array. 

[00059] FIGs. 3 A-F show some exemplary incident wavefronts and corresponding output 
irradiance distributions in the focal plane. 

[00060] FIGs. 4A-F show some additional exemplary incident wavefronts and corresponding 

output irradiance distributions in the focal plane. 

[00061] FIG. 5 shows simulated CCD data. 

[00062] FIGs. 6A-B illustrate center of mass calculations. 

[00063] FIG. 7 plots error in the center of mass calculation. 

[00064] FIGs. 8-25 illustrate other aspects of the invention. 

DETAILED DESCRIPTION 
[00065] OPTIMIZED METHODS FOR FOCAL SPOT LOCATION USING CENTER OF 
MASS ALGORITHMS. 

[00066] We want to know the location of the focal spot and we use the center of mass technique 
to approximate this location. The question is: how well does the center of mass approximate the 
focal spot location? To resolve this issue we begin with a special case, and then move to the 
situation encountered in the laboratory. 

[00067] Consider the special case of a spatial irradiance function g(x) being symmetric about the 
peak in some local region. In other words, the function has even parity about the peak xO in some 
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local region which is smaller than the integration domain for the center of mass. The requirement 
for even parity can be written as: 

g(x-Xo)~g{x+x 0 ). (25) 
Since g(x) is even, the antiderivative of this function is odd on this same interval. Defining the 
antiderivative as: 

lg(x)dx-G(x), (26) 
the parity condition becomes: 

G(x-x 0 )--G(x+x 0 ). (27) 

[00068] Now we consider the center of mass computation. Defining the center of mass on some 
domain [a, b] as: 

\xg(x)dx 

x=\ (28) 

\g(x)dx 

we now have the most general formula. But because of the parity of g(x), we are motivated to 
consider a symmetric domain centered on the peak. This implies: 

x 0 -a=-(x 0 -b). (29) 
[00069] If we consider p to be radius we can define: 

p=x 0 -a (30) 
and the center of mass in equation 28 now becomes: 
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2*.<KxJ%r,- \G(x)dx 



X = 



(31) 



[00070] Since G(x) is odd over this domain, the integral vanishes. And using equation 27 the 
center of mass can be reduced to: 



[00071] So when g(x) is symmetric about the peak and the domain is centered on the peak, the 
center of mass is exactly the focal spot location. 

[00072] Due to the discreet nature of the CCD array, the pixel addresses are integers and this 
precludes one from selecting integration boundaries that are centered exactly about the peak. This 
leads to the consideration of how the center of mass shifts away from the peak. The approach 
here is to consider the problem as having two parts: a piece with exact symmetry and the piece 
which shifts the center of mass away from the focal spot location. 

[00073] For the arbitrary domain [a, b], no longer centered about x 0 , the center, c, of the domain 
is: 



[00074] This leads to the definition of an asymmetry parameter x which measures how far the 
peak is from the domain center: 




(32) 



a+b 



2 



(33) 



T = X 0 -C. 



(34) 



[00075] The next step is to define the radius p of the symmetric piece: 
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p= Min(x 0 -a,b- x Q ). (35) 

[00076] The center of mass can now be written in terms of the symmetric and asymmetric 
component as: 

Xq+P x 0 +p-2t 



jxg(x)dx+ jxg(x)dx 

~_P Xo+p 

o+P x 0 +p-2r 

\g(x)dx+ \g{x)dx 

x 0 -p x 0 +p 



A " x 0 +p x 0 +p-2T ' (30) 



which reduces to: 



x 0 +p-2r 



(x 0 - p)G(x 0 + p)+ (x 0 + p- 2t)G(x q + p-2r)- j G(x)dx 

x= x ^- p 



(37) 



G(x 0 + p)+(x 0 + p-2r) 

[00077] As x — »0 the integration boundaries become symmetric about the peak and the center of 
mass approaches the peak, i.e. x Xo. 

[00078] The error 8 = x - x 0 in the center of mass result can be expressed in terms of the 
asymmetry parameter x which is a key to a subsequent correction scheme. For the case of a 
sinc2x function, this error is given as: 

lm»- ln(p- 2r)+ ci{l(p- 2t))- ci(2p) 
2 3Si(2p)-sh(p- 2r)) + Sm{p - 2T > - 3^- 

where ci(x) is the cosine integral function: 

f»cosr 

ci(x) =-\ x —dt , (39) 
[00079] and Si(x) is the sine integral function: 
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&■(*) = J— Tft, (40) 
o ' 

[00080] We now present three correction strategies to improve the center of mass computation 
and recover a more accurate focal spot location. 
[00081] Integration Domain Boundary Method. 

[00082] One method for implementing the center of mass to exploit the results of the analysis 
above locates the brightest pixel under a lenslet and then marches out in the ±x and ±y directions 
to find the four minima. These four pixels define the boundary of a rectangular integration 
domain. The smallest of the four minima is used as a threshold and this value is subtracted from 
all pixels in the integration domain. 

[00083] FIG. 5 shows simulated CCD data for illustrating how an irradiance distribution is 
represented discretely. Here the incident wavefront was isoplanatic and with a small tilt 
component. The rectangular box defined by the white lines represent the integration region. Here 
the error is pixels. 6 T = (0.013, 0.002) 

[00084] FIGs. 6A-B shows that the center of mass calculation does not match the focal spot 
location in the general case (FIG. 6A) due to an asymmetry. The center of mass can be thought 
of as having two components: a symmetric piece (light gray) and the asymmetric component 
(dark gray). FIG. 6B shows that removing the pedestal greatly reduces the shift caused by the 
center of mass computation by reducing the contribution of the asymmetric piece. The 
asymmetric piece is the contaminant that pulls the center of mass value away from the peak 
location. Successful center of mass strategies will minimize the effect of this piece. 
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[00085] FIGs. 6A-B illustrate the wisdom of the Integration Domain Boundary method. First of 
all, it endeavors to symmetrically bound the peak, which minimizes the error in the center of 
mass computation. This also gobbles up more pixels than a thresholding method, which provides 
more information. Now the asymmetric contaminant is also minimized near the minimum. The 
total area in the darker region is what causes the error in the computation. Near the minimum, 
this piece has the smallest height. And after removing the pedestal (the dark count from the 
camera), the piece is even smaller. The width of the piece is now less than the width of a pixel. 
So this simple scheme exploits the lessons from the analysis above all without having to provide 
a threshold a priori. 

[00086] There is a drawback to this method in that as the distance between the asymmetric piece 
and the peak increases, the lever arm also increases. In other words, a unit area added to the 
center of mass is more harmful at larger distances. However it appears that in the cases of square 
lenslets with well-defined sinc2x sinc2y focal spots, the lever arm effect is more than 
compensated for by the height minimization at the minimum. This seems to suggest against using 
the second minimum as the lever arm doubles. 

[00087] Some basic numerical studies were on simulated camera files created using the 
propagation equations defined above. These data did not include either noise or cross-talk from 
neighboring lenslets. The advantage of simulated data is that one has exact knowledge of the 
peak and can exactly quantify the error in the center of mass calculation. For comparison 
purposes, a classic center of mass with a threshold of 12% was used. The average error 5 r = (5,, 
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Sy) for each method was then computed for 4096 samples all where the peaks were within one 
pixel of the projection of the lenslet center. The standard deviation of the errors are: 
crj = {0.015, 0.015 } pixels for the pure center of mass method, and 
crj = {0.010, 0.01 1 } pixels for the integration domain boundary method. 
[00088] Modified Integration Domain Boundary Method. 

[00089] The above Integration Domain Boundary method is an excellent starting point for 
improving the center of mass calculation. In the cases where the form of the irradiance 
distribution is known, an additional improvement is possible. Since the shape of the distribution 
is known, one can look at the difference between the center of mass and the asymmetry 
parameter and make a quick correction based on a two-dimensional version of equation (38). In 
the sample environment tested, the additional correction was able to remove 90% of the error, 
reducing the error by an order of magnitude. 

[00090] FIG. 7 shows the error in the center of mass computation as a function of peak offset 
from the center of the pixel. The axes describe the distance the peak was offset from the center of 
the pixel. The color of the block reveals the magnitude of the error. As expected, there is no 
error when the peak is centered in the pixel. The maximum error in this example was 0.03 pixels. 
The effect of the discrete change in integration boundaries manifests distinctly as an abrupt 
change in the error. 

[00091] The error variations in FIG. 7 were studied in the x and y components. The next step 
was to plot the error components as a function of the asymmetry parameter. A surprisingly 
simple variation was observed and is shown in FIG. 8. FIG. 8 plots the center of mass error as a 
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function of the asymmetry parameter. Here the graph shows the x component of the error plotted 
in FIG. 7. Note that the functional relation between error and the asymmetry parameter is linear 
with two distinct branches. The errors are essentially linear with respect to the asymmetry 
parameter x with two distinct branches. This implies that a simple linear correction can be used to 
push the center of mass closer to the focal spot location. For these data the results of the linear 
fits were: 

positive branch error = (0.043756 ± 0.000004) x + (0.00000 ± 0.00002) 
negative branch error = (-0.03055 1 5 ± 0.000004) x + (0.00000 ± 0.00002) 
[00092] The beauty of this method is that the slopes for the positive and negative branches are 
functions of the wavefront sensor. As such, they only need be determined one time. The method 
is to compute a center of mass and an asymmetry parameter. Then determine which branch of the 
correction is appropriate, and apply a simple linear boost to nudge the center of mass closer to the 
peak location. The nudge is applied separately to the x and y components of the center of mass 
prediction. 

[00093] FIG. 8 allows one to see graphically how precise this correction can be. The vertical 
spread in each cluster of points is approximately 0.002 pixels. The linear correction basically cuts 
this in half, leaving us with an error of approximately 0.001 pixels. Attempts to reduce the error 
further by handling the problem in two dimensions failed. 

[00094] For example, if the asymmetry parameter was at the maximum of x = 0.3 pixels, the 
correction is: 

(0.043756 ± 0.000004)0.3 + (0.00000 ± 0.00002 ) = -0.01313 ± 0.00002 pixels. 
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which is added to the computed center of mass. This will bring the computation to within 0.001 
pixels of the true focal spot location. 

[00095] Another strategy has been developed in collaboration with J. Roller. This is an iterative 
scheme. The gist of the method is that one computes the center of mass and develops an 
approximate focal spot location, (x 0 , yd). This location should easily be within a pixel of the true 
focal spot location. A new integration domain is now defined, centered on the center of mass 
value. The size of the domain is given by the Tsvetkov criteria: use all the data within the first 
minima. We know from equation 9 that the half-width is p = fk/d. So the integration domain is 
[x 0 - p, x 0 + p] and the integration range is [y 0 - p, yo + pL Using these new boundaries, a new 
center of mass is computed, (x\ 9 y\) 9 along with new boundaries. The process continues until the 
distance between peak locations on consecutive iterations is smaller than some arbitrary 
acceptance parameter. 

[00096] The problem is that now the integration boundaries are no longer integer pixel values, 
but they lie within the pixels. So the heart of this scheme is to find a way to apportion the light 
within a pixel and maintain the speed advantage of the center of mass. While the exact 
computation of this apportionment can be messy, the exact solution can be approximated quite 
well with a polynomial. 

[00097] Consider the cases of sine 2 * in one dimension. The parameter cp describes where the 
minimum falls within the pixel. For example, if the minimum is in the middle of the pixel, <p - 
0.5; if the minimum is three quarters of the way across the pixel then <p - 0.75. The function co((p) 
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is weighting function varying from zero to one which describes how much of the light to include 
from the pixel. In this case the weighting function is given exactly as: 
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(41) 



[00098] The constant a represents: 

a = Si(2K)~ 1.4181515761326284502457801622997494291424533492950 .... (42) 
[00099] Notice that this is for pixels on the right-hand side. One can exploit the even parity of 
the sine 2 * function to get the values for the left-hand side. 

[000100] The problem is of course that a naive implementation of Si(x) will slow down the 
algorithm tremendously. We can exploit the fact that over the domain interest, 9 e [0,], the 
function has a rather simple behavior. FIG. 9 shows the weighting function as a solid curve. The 
points are the results of polynomial least squares fit through order nine. The agreement is quite 
good. 

[000101] We can write the approximate solution Q((p) as: 



(43) 



[000102] wherein n represents the degree of fit. In FIG. 9, n = 9. Notice that this formulation 
avoids the ambiguous quantity 0°. Of course a direct implementation of equation 43 would be 
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ponderously slow and erode the speed advantage of the center mass technique. We show two 
different ways to encode these results in C++: 



double 


Ctl0] ; 


//LSF fit coefficients 


c[0] = 


783 .7320568947441; 


//constant term 


c[l] = 


-3448 .7204133862506; 


//linear term 


c[2] = 


6172.40085200703; 


//coefficient for x^2 


c[3]' = 


-5739.328837072939; 




C[4] = 


2945.8812707469388; 




c[5] = 


-833.7602125492338; 




c[6] = 


130.30391297857992; 




c[7] = 


-9.792719111240226; 




c[8] = 


0 .2907165577816425; 


//coefficient for 



double omega (double * c, double phi, int n) //weighting function 



// c LSF fit coefficients for equation 41 

// phi pixel fraction 

// n degree of polynomial least squares fit 

// lsf weighting factor to use (omega [phil] ) 

{ 

double lsf = 0 //LSF prediction 

for (int i = n; i > 0; i--) //loop through orders 

lsf = phi * (lsf + c[i]); //build up answer 
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lsf += c[0]; //add in constant 

term 

return lsf; //return the answer 

} 

[000103] This is an elegant implementation that easily accommodates fits of arbitrary order 
through adjusting the parameter n. However, a somewhat faster method is 

double f (double * c, double phi) //rapid polynomial 

{ 

double lsf = c[0] + phi * (c[l] + phi * (c [2] + phi * (c[3] 

+ phi * (c[4] + phi * (c[5] + phi * (c[6] 
+ phi * (c[7] + phi * (c[8] + phi * 

(c[9]))))))))) ; 
return lsf; 

} 

[000104] The first method took 16% longer to execute than the second method. A 
straightforward implementation of equation 43 took 13 times longer than the second method. 
[000105] We have shown formally the well-known result that the shift in the focal spot location 
is proportional to the average wavefront slope across a lenslet. The result is exact in the 
isoplanatic and small angle limits. We then showed the relationship between the focal spot 
location and the center of mass computation. The accuracy of the assumption that the center of 
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mass is the focal spot location depends directly upon how symmetrically the integration domain 
bounds the peak. When the domain is exactly symmetric, the center of mass exactly represents 
the peak. Of course with real data, the integration boundaries are not centered upon the peaks. 
We closed with three strategies to improve the accuracy of the center of mass computation and 
preserve the time advantage of the method. 

[000106] PRECISION MEASUREMENTS OF A LENSLET ARRAY. 
[000107] As discussed above, a critical step in data processing involves fixing the focal spot 
location and a common method used for quickness is the center of mass. But naive applications 
of this method have a precision of 0.01 pixels, which forms an inherent precision boundary for 
all computations. 

[000108] While the error in the center of mass technique is highly systematic, its manifestation 
in an experimental measurement is typically, to first order stochastic. The challenge then 
becomes how to remove these stochastic errors. The first tool in such situation is the least 
squares fit which is a simple way to remove stochastic errors. 

[000109] Here we describe how a least squares fit to used to precisely measure the size and 
location of a lens-let array. Here we assume a rectangular lenslet and slide the grid around to 
find the best location assuming the focal spot locations are under the center of each lenslet. 
[000110] The least squares fit. 

[000111] The first decision is which parameters to extract from the least squares fit. Although 
we could assume a square lenslet configuration, we used a rectangular configuration. Of course 
we expect to find a square, but by solving for the two widths (dx, dy) independently, we are able 
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to check the quality of the result by seeing if the aspect ratio is unity. Also, this formulation 
allows for the more general case of rectangular lenslets. The other free parameter is the location 
of the lenslet origin, O = (xo, yo). 

[000112] The fitting strategy is elementary. Consider two rows of points, the points 
corresponding to the focal spot locations under each lenslet. We would want to find the height of 
a line which best bisects the two rows. 

[000113] Here r is the row index and c is the column index andjy is the abscissa of the focal 
spot location. The summation is across all columns and involves and arbitrary row r and the next 
row, r + 1 . Clearly, if all the points rested upon two horizontal lines, the merit function x 2 would 
be zero when h is the bisector. 

[000114] This method is easily extended to the other dimension by finding the line which best 
bisects each column. To formulate the full problem in two dimensions, we need to define a few 
quantities. We already have the origin vector 0, 



o= 



yo. 



(45) 



and the lenslet width vector d, 



d = 



dx 
dy 



(46) 

[000115] We need to define a position vector i, which specifies the row and column the lenslet: 
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(47) 



[0001 16] Finally we note that the focal spot location for the lenslet in column c and row r of the 
lenslet array is: 



(48) 



so we are loading the x andy locations into a vector for efficient computation. We note that the k 
vector describes the number of rows and columns in the lenslet array that are projected onto the 
CCD array: 



[000117] Finally we need a collection of basis vectors. In two dimensions, it is a spinor: 



(49) 



(50) 



[000118] This structure points to the neighboring object. Although this is a two-dimensional 
problem, this methodology can be trivially extended to an arbitrary number of dimensions and it 
reduces the complexity of computer code. 

[000119] If we use as a spacetime index to sum over the dimensions, the full merit function is 

2 k -°„ 



„=1 /»=(!,!) V ZV 



(51) 
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[000120] The problem is solved in the canonical fashion. All first derivatives of the fitting 
parameters are forced to zero simultaneously by solving a set of linear equations. That is, the 
four simultaneous requirements that: 

doX 1 = 0, (52) 



and 



^ 2 = o, 



(53) 



[000121] generate a series of matrix questions. The solution to these equations follows. 
[000122] First, we define some intermediate variables. There are two summation variables, n 
and a. They are given by: 



k-a u 



n » = X 1 > 



(54) 



and 



k-<r u 



[000123] If every lenslet registers a focal spot, then these quantities are (CRC Standard 
Mathematical Tables, p 54): 



(55) 
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a = 
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(56) 



k, 



(57) 



and 
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(58) 



[000124] Finally, we define a determinant, A: 

[000125] We can now begin writing out the solutions. The lenslet origin is: 

k-a„ 



(59) 



4,= 



2A 



(60) 



[000126] The lenslet size is: 

k-o..( 



*-<r„ 



2A 



(61) 



[000127] We now turn our attention to the computation of statistical uncertainties. The first 
quantity to compute is the sample standard deviation (Bevington, p 19): 

2 k ~o M 



2 p=(l,l) 

S = 



n M -2 



(62) 



[000128] This allows us to compute the uncertainties in the size of the lenslets and the position 
of the array. The uncertainty in the lenslet origin is: 

\2\ 2 

(63) 



^2 _ 2 
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A 
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and the uncertainty in the lenslet width is: 




(64) 



[000129] A sample lenslet and camera were chosen from the Columbus experimental 
configuration. This sensor uses a 60 x 60 array of square lenslets. Each lenslet is 280 |i across 
with a spherical profile and a focal length of 28 mm. The array is etched into fused silica and the 
lenslets have a 100% fill factor. The camera used was an SMD model 1M15, with a 14 |i pixels 
in a 1024 x 1024 grid. The output is digitized to 12 bits. The lenslet array was placed 
approximately 36 mm (1.3 focal lengths) from the CCD array. 

[000130] By design, the lenslet array is much larger that the CCD array, so a subset of the 
lenslets are used. Here we used a 50 x 50 subset of the lenslets. For these 2500, the array, was 
found to be: 

size = 20.0005 ± 0.0008 x 20.0025 ± 0.0008 pixels, 
and location of the origin of the array was found at the pixel address of: 

location = (16.92 ± 0.02, 15.58 ± 0.02) pixels. 
[000131] The aspect ratio of the lenslet is then: 



implying that the lenslets are very close to being square. 

[000132] Of course this measurement does not simply measure the regularity of the lenslet 
array, it also measures the regularity of the CCD array. But the basic technique demonstrates 



aspect = 



20.0025+0.0008 
20.0005 ±0.0008 



= 1.000101 0.00006 



(65) 
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clearly that an instruments with an inherent low precision limit can be an effective tool to make 
precise measurements. 

[000133] The data above are part of a larger set of data collected by T.D. Raymond. An optical 
fiber was set 750 mm from a collimating lens and a shear plate was used to check collimation. 
These data are the baseline file, xOyO.bvf. Then a differential micrometer was used to translate 
the fiber in a plane perpendicular to the direction of propagation. This method introduced up to 
400 nr of tip and tilt in 50 nr increments. The data were analyzed for each position and a subset 
of the results are shown in table 1 below. 

[000134] A close look at the data reveals a very slight asymmetry, at the fringe of the 
experimental uncertainty. Most formally this number is ratio of the aspect ratio of the lenslet to 
the aspect ratio of the CCD array. The natural question becomes is the asymmetry in the lenslet 
array or the CCD array? 



Table 1 



File 


origin 


size 


Aspect 


xOyO 


16.92 ±0.02, 15.58 ±0.02 


20.0005 ± 0.0008 x 20.0025 ± 0.0008 


1.00010 ±0.00006 


x0y50 


16.79 ±0.02, 15.58 ±0.02 


20.0006 ± 0.0008 x 20.0026 ± 0.0008 


1.00010 ±0.00006 


xOylOO 


16.67 ±0.02, 15.58 ±0.02 


20.0005 ± 0.0008 x 20.0027 ± 0.0008 


1.00011 ±0.00006 


x0yl50 


16.55 ±0.02, 15.58 ±0.02 


20.0005 ± 0.0008 x 20.0026 ± 0.0008 


1.00011 ±0.00006 


x0y200 


16.41 ±0.02, 15.58 ±0.02 


20.0007 ± 0.0008 x 20.0026 ± 0.0008 


1.00010 ±0.00006 


x0y250 


16.29 ±0.02, 15.58 ±0.02 


20.0006 ± 0.0008 x 20.0027 ± 0.0008 


1.00010 ±0.00006 


x0y300 


16.16 ±0.02, 15.58 ±0.02 


20.0006 ± 0.0008 x 20.0027 ± 0.0008 


1.00010 ±0.00006 


x50y0 


16.92 ±0.02, 15.70 ±0.02 


20.0007 ± 0.0008 x 20.0030 ± 0.0008 


1.00011 ±0.00006 
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JT lie 


ungin 


size 


Aspect 


x50y50 


16.80 ±0.02, 15.70 ±0.02 


20.0005 ± 0.0008 x 20.0029 ± 0.0008 


1.00012 ±0.00006 


x50yl00 


16.67 ±0.02, 15.70 ±0.02 


20.0005 ± 0.0008 x 20.0029 ± 0.0008 


1.00012 ±0.00006 


x50yl50 


16.55 ±0.02, 15.70 ±0.02 


20.0006 ± 0.0008 x 20.0029 ± 0.0008 


1.00011 ±0.00006 


x50y200 


16.41 ±0.02, 15.70 ±0.02 


20.0008 ± 0.0008 x 20.0030 ± 0.0008 


1.00011 ±0.00006 


x50y250 


16.29 ±0.02, 15.70 ±0.02 


20.0007 ± 0.0008 x 20.0029 ± 0.0008 


1.00011 ±0.00006 


x50y300 


16.17 ±0.02, 15.70 ±0.02 


20.0006 ± 0.0008 x 20.0028 ± 0.0008 


1.00011 ±0.00006 


xlOOyO 


16.92 ±0.02, 15.82 ±0.02 


20.0007 ± 0.0008 x 20.0032 ± 0.0008 


1.00012 ±0.00006 


xl00y50 


16.80 ±0.02, 15.82 ±0.02 


20.0005 ± 0.0008 x 20.0030 ± 0.0008 


1.00013 ±0.00006 


xlOOylOO 


16.67 ±0.02, 15.82 ±0.02 


20.0006 ± 0.0008 x 20.0032 ± 0.0008 


1.00013 ±0.00006 


xl00yl50 


16.55 ±0.02, 15.82 ±0.02 


20.0006 ± 0.0008 x 20.0030 ± 0.0008 


1.00012 ±0.00006 


xl00y200 


16.42 ±0.02, 15.82 ±0.02 


20.0006 ± 0.0008 x 20.0030 ± 0.0008 


1.00012 ±0.00006 


xl00y250 


16.29 ±0.02, 15.82 ±0.02 


20.0008 ± 0.0008 x 20.0032 ± 0.0008 


1.00012 ±0.00006 


xl00y300 


16.16 ±0.02, 15.82 ±0.02 


20.0007 ± 0.0008 x 20.0033 ± 0.0008 


1.00013 ±0.00006 


xl50y0 


16.92 ±0.02, 15.95 ±0.02 


20.0006 ± 0.0009 x 20.0028 ± 0.0009 


1.00011 ±0.00006 


xl50y50 


16.79 ±0.02, 15.96 ±0.02 


20.0006 ± 0.0009 x 20.0029 ± 0.0009 


1.00012 ±0.00006 


xlSOylOO 


16.67 ±0.02, 15.95 ±0.02 


20.0006 ± 0.0008 x 20.0031 ± 0.0008 


1.00012 ±0.00006 


xl50yl50 


16.55 ±0.02, 15.96 ±0.02 


20.0006 ± 0.0008 x 20.0028 ± 0.0008 


1.00011 ±0.00006 


xl50y200 


16.42 ±0.02, 15.96 ±0.02 


20.0006 ± 0.0008 x 20.0027 ± 0.0008 


1.00011 ±0.00006 


xl50y250 


16.29 ±0.02, 15.96 ±0.02 


20.0007 ± 0.0008 x 20.0028 ± 0.0008 


1.00010 ±0.00006 


xl50y300 


16.18 ±0.02, 15.95 ±0.02 


20.0004 ± 0.0008 x 20.0028 ± 0.0008 


1.00012 ±0.00006 


x200y0 


16.92 ±0.02, 16.10 ±0.02 


20.0007 ± 0.0009 x 20.0022 ± 0.0009 


1.00007 ±0.00006 


x200y50 


16.80 ±0.02, 16.10 ±0.02 


20.0006 ± 0.0008 x 20.0022 ± 0.0008 


1.00008 ±0.00006 


x200yl00 


16.68 ±0.02, 16.10 ±0.02 


20.0005 ± 0.0008 x 20.0022 ± 0.0008 


1.00009 ±0.00006 
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rue 


origin 


size 


Aspect 


x200yl50 


16.56 ±0.02, 16.10 ±0.02 


20.0005 ± 0.0008 x 20.0022 ± 0.0008 


1.00008 ±0.00006 


x200y200 


16.42 ±0.02, 16.10 ±0.02 


20.0007 ± 0.0008 x 20.0022 ± 0.0008 


1.00007 ±0.00006 


x200y250 


16.29 ±0.02, 16.10 ±0.02 


20.0007 ± 0.0008 x 20.0022 ± 0.0008 


1.00007 ±0.00006 


x200y300 


16.18 ±0.02, 16.10 ±0.02 


20.0005 ± 0.0008 x 20.0020 ± 0.0008 


1.00007 ±0.00006 


x250y0 


16.93 ±0.02, 16,23 ±0.02 


20.0007 ± 0.0008 x 20.0018 ± 0.0008 


1.00006 ±0.00006 


x250y50 


16.80 ±0.02, 16.23 ±0.02 


20.0006 ± 0.0008 x 20.0017 ± 0.0008 


1.00005 ±0.00006 


x250yl00 


16.67 ±0.02, 16.22 ±0.02 


20.0007 ± 0.0008 x 20.0019 ± 0.0008 


1.00006 ±0.00005 


x250yl50 


16.54 ±0.02, 16.22 ±0.02 


20.0007 ± 0.0008 x 20.0020 ± 0.0008 


1.00006 ±0.00005 


x250y200 


16.41 ±0.02, 16.23 ±0.02 


20.0008 ± 0.0008 x 20.0019 ± 0.0008 


1.00005 ±0.00005 


x250y250 


16.30 ±0.02, 16.23 ±0.02 


20.0006 ± 0.0008 x 20.0018 ± 0,0008 


1.00006 ±0.00005 


x250y300 


16.17 ±0.02, 16.24 ±0.02 


20.0006 ± 0.0008 x 20.0018 ± 0.0008 


1.00006 ±0.00005 


x300y0 


16.93 ±0.02, 16.35 ±0.02 


20.0007 ± 0.0008 x 20.0021 ± 0.0008 


1.00007 ±0.00006 


x300y50 


16.80 ±0.02, 16.35 ±0.02 


20.0006 ± 0.0008 x 20.0021 ± 0.0008 


1.00008 ±0.00006 


x300y!00 


16.69 ±0.02, 16.35 ±0.02 


20.0004 ± 0.0008 x 20.0020 ± 0.0008 


1.00008 ±0.00006 


x300yl50 


16.56 ± 0.02, 16.35 ± 0.02 


20.0005 ± 0.0008 x 20.0020 ± 0.0008 


1.00007 ±0.00006 


x300y200 


16.42 ±0.02, 16.35 ±0.02 


20.0005 ± 0.0008 x 20.0022 ± 0.0008 


1.00007 ±0.00005 


x300y250 


16.30 ±0.02, 16.35 ±0.02 


20.0007 ± 0.0008 x 20.0021 ± 0.0008 


1.00007 ±0.00005 


x300y300 


16.17 ±0.02, 16.35 ±0.02 


20.0007 ± 0.0008 x 20.0022 ± 0.0008 


1.00008 ±0.00005 



[000135] A differential micrometer was used to introduce tip and tilt in 50 ur increments. The 
computations in equations 60-64 were then performed on the new data set. 
[000136] WAVEFRONT RECONSTRUCTION FOR THE SHACK-HARTMANN 
WAVEFRONT SENSOR. 
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[000137] FIGs. 2A-F show a Shack-Hartmann sensor tracking a series of focal spots. This 
idealization shows the process of measuring a spherical wave on a sensor with just 16 lenslets. 
The first step, as represented by the FIGs. 2A, 2C and 2E, is to measure a plane wave and 
measure a series of focal spot locations which are used as a reference. The next step, as depicted 
in FIGs. 2B, 2D and 2F, is to introduce the wavefront of interest and determine the focal spot 
shifts. The critical information at this stage of the analysis is the focal spot shifts for every 
lenslet. 

[000138] As discussed above, one uses the results of myriad averaged slope measurements to 
reconstruct the wavefront incident upon the lenslet array. It is common to hear people describe 
reconstruction as an integration process. This is malapropos and ironic. The irony is that Nature 
has already done the integration for us; as shown below the measured quantity is already an 
integrated form of the incident wavefront. 

[000139] Given a wavefront y (xl, x2) incident upon a lenslet, what exactly does the sensor 
measure in the isoplanatic and small angle limits? It measures the average of the wavefront 
slopes ^ over each lenslet. So for an integration domain of xi e [a, b], x 2 e [c, d], the average 
of the [i derivative of \|/(xl, x2) is given by: 



where the operator dp. denotes differentiation along a line parallel to the |li axis. Defining the 
antiderivative of ^(x 1 ^ 2 ) as: 



d b 




c a 



(66) 



(d-c)(b-a) 
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M*i'*2)= lw(x {9 x 2 )dx M , (67) 
the average of the slope is written as: 

= A- l (y, v (a 9 c)- y v (a 9 d)- y v (b,c)+ y v (b,d)) (68) 
where the domain area has be recast as A = (d - c)(b - a) and v represents the orthogonal 
direction. The message of equation 68 is quite clear: the measurement of the focal spot shift tells 
us about differences in the antiderivatives evaluated at the lenslet vertices. 
[000140] At this juncture, the reconstruction will typically follow one of two paths: either 
modal or zonal reconstruction. An example of the two different types of outputs is shown in 
FIGs. 10-13B. 

[000141] Zonal methods output a value \|/ at each zone. In the Shack-Hartmann case, the 
wavefront values are computed at the vertices. In the Southwell case, the wavefront value is 
computed for the center of the lenslet. The modal reconstructor output the amplitudes of some 
arbitrary set of functions used to describe the wavefront. A practical set of functions are the 
Taylor monomials. With these amplitudes, it is trivial to compute any other rational function set. 
Of course, one can not easily convert the Taylor amplitudes into amplitudes for a irrational 
function set, such as one containing Gaussians. 

[000142] Consider a fit to order d using a sensor with an § x r| lenslet array. The output from 
the zonal reconstruction will be + l) 2 (r| + l) 2 wavefront values for each of the vertices, 
reducing the wavefront to a discrete data set. The modal reconstructor will yield (d + l)(d + 2)/2 
amplitudes, and the wavefront is now in a continuous representation . As an example, consider a 
6 th order fit with a 50 x 50 lenslet array. The zonal fit will yield 2601 numbers; the modal fit will 
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produce 28 numbers. So the information content of the zonal reconstruction is typically richer. 



One may expect to pay a computation burden for this extra information, but this is not 
necessarily the case as we shall see. 

[000143] FIGs. 10-13B show the reconstruction process schematically. In this example, an 
eigenstate wavefront of the Zemike polynomials, Z42, has been imaged onto a 48 x 48 lenslet 
array. The focal spot shifts are computed and passed to either the zonal or the modal 
reconstructor. A sample output from each reconstructor is shown. The zonal reconstructor 
provides a value for the wavefront at every lenslet vertex (FIG. 13 A); the modal reconstructor 
provides an amplitude for every function in the basis set (FIG. 13B). 

[000144] A wavefront is imaged onto the lenslet array creating a pattern of focal spots on the 
CCD array. The software then computes how far the focal spots have moved from the reference 
positions and creates a series of offset measurements, 5. The focal spot shift is related to the 
lenslet focal length / and the average wavefront slope m by: 



6-f 



m 



(69) 



when the CCD is in the focal plane of the lenslet array. In practice, the radical in the 



denominator is so close to unity that it is ignored leading to: 



b=fin. 



(70) 



[000145] The solution to the full form is: 




(71) 



and the approximation is: 
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8 

m=j. (72) 

[000146] For example, an array with a focal length of 20 mm is imaged onto a CCD array with 
24|x pixels. A shift of two pixels is measured. The wavefront slope then is 1 .2 mr using the 
approximation in equation 72, and is 1.200000864 mr using the exact form in equation 71. So 
for practical purposes, the approximation is adequate for lenslet arrays with long focal lengths. 
[000147] The measurement of the focal spot shift 5 then is a measurement of the average slope 
over the lenslet. In the isoplanatic limit the average slope can be thought of as a tilt. The effect 
then of a measurement is to reduce the wavefront into a series of planar patches of the size of 
each lenslet, with a tilt in each direction corresponding to the focal spot shift in that direction. 
[000148] To delve deeper into the details of the calculation, we need to introduce some 
notation. There are a series of N offset measurements where N= fy\ is the number of lenslets. 
Eventually, the linear algebra of zonal reconstruction is going to force the lenslets to be ordered 
in a linear fashion. In the interest of expediency, we will use the same ordering for both the 
zonal and the modal cases. 

[000149] The ordering of the lenslets is completely arbitrary: it can be row major, column 
major or even random. The position vector p describes the linear ordering of the addresses. A 
sample p vector is: 

/> = ((!, 1),(2,1), (3,1)...) (73) 
where the address pairs (w, v) denotes rows and columns. On an £ x r| array, the address (u, v) 
would be at position w(£ - 1) + v in the p vector. The reverse look-up procedure is a bit more 
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involved. It is simple to compute, but in practice, it is convenient to construct a lookup able q 
which is really an £ x r| matrix. A sample q matrix is shown below. 



1 

7 + 1 



2 

7 + 2 



L(£- 1)7+1 (£-1)7 + 2 



27 



<?7 



(74) 



[000150] Equation 68 has already paved the road for modal reconstruction. Before beginning 
with the details of reconstruction, it is helpful to define the box operator for an arbitrary 
function over the domain x\ e [a, b], x 2 e [c, d]. This operator combines the values of a function 
at the domain corners in the following way: 

V(x) = y(a, c) - y(a, d) - y(b, c) +\|/(6, d). (75) 
[000151] Having seen the box operator in an explicit form, we sill now generalize its 
formulation. Consider the rectangular domain show in FIG. 22. FIG. 22 shows the vertices of a 
lenslet. Here the vertex points are labeled with the Greek letters a through 5. Each of these 
letters is really an address pair. The indexing of the vertices is tied to the indexing for the lenslet. 
If this is lenslet (u, v), then a is vertex (u-\ s v-lj; p is u, v-7), y is (u, v), and 5 is (u-l, v). The 
vertices are labeled by a, P, y, and 6 where each vertex has an address (x u x 2 ). These vertices can 
be arranged in a vertex vector V as: 



(76) 



[000152] When a function operates on this vector, it returns a vector of functions. Explicitly, 
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(77) 



[000153] (To save space, subsequent vectors will be written in the transpose form.) This leads 
to a generalization of operators which act upon the vertices. Think of an interaction matrix T 
which describes the scalar interaction of the vertices. For the box operator: 



-10 0 0 
0 10 0 
0 0 10 
0 0 0 1 



(78) 



and we can now write 



v|/(x) = r fl v|/(F). (79) 
[0001541 With this box operator, the average of the (i slope over the lenslet at address /?, is 
expressed as: 



(80) 



[000155] In this case (equation) is the predicted value. But the actual measurement is of 
(equation). To find the wavefront which best describes the measurements, one simply needs to 
minimize the merit function: 



fj=l 1=1 



(81) 
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where a,- represents an error estimate. However, since these errors do not depend upon the 
amplitudes, they can be treated as weighting factors and are withheld from explicit evaluation. 
[000156] All that remains now is to choose a set of basis functions. This set should be 
complete. Orthogonal is sets are attractive, but often involve unneeded overhead. The set we 
will choose is the basis set for all rational functions: the Taylor monomials. These monomials 
are given by: 

r- jJ {x x ,x 2 )=x[- j x J 2 (82) 
and the wavefront described as a c/th order expansion would be written as: 

*W=ZI« w (83) 

i-O y-0 '"J J 1 2 

[000157] We specify this set because it contains all the information needed to construct the 
functions of all other rational basis sets, such as the polynomials of Zernike or Legendre. It is 
easiest to work in the Taylor basis and then use an affine transformation on the Taylor amplitudes 
to compute the amplitudes in the other bases. These transformations matrices are written in 
terms of integers, so there is no loss of precision. 

[000158] Also, the Taylor basis far easier to work with. Consider the case of the 
antiderivatives. They can be written to arbitrary order as: 

t JJ (x x ,x 2 )=j-j^-x\-Jxi, (84) 

and 

n- JJ ( x ^h-}]-< J ^ +l (85) 
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[000159] By comparison, writing out the two antiderivatives for each of the 66 Zermike 
polynomials through 10 th order would fill several pages. 

[000160] Anyway, the following analysis does not depend upon the basis set. The 
minimization of the merit function proceeds in the canonical way: the derivatives with respect to 
each amplitude are set to zero, 



[000161] This will generate a series of n = (d + \){d +2)/2 linear equations to solve. Since this 
is a linear least squares fit, we expect to find a direct solution. All that remains are some linear 
algebra details unique to the basis set chosen. In the interest of brevity, we will focus to the issue 
of zonal reconstruction. 

[000162] We now consider the problem of zonal reconstruction for the Shack-Hartmann sensor. 
Here too we start with equation 68 and derive the merit function A (which is %2 without the a, 
terms): 



[000163] There is a problem with the antiderivative terms. They are quite unwelcome for a 
zonal fit because it would require discrete differentiation to recover to wavefront values. It is 
better to recast A to explicitly compute the wavefront values. So we turn our attention to the task 
of eliminating the antiderivatives. We begin with a close examination of equation (66). Note 
that in one dimension, the average value of the derivative of a function g(x) over the linear 
domain xl G [a, b] is trivial to compute: 




0. 



(86) 




(87) 
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0 \S\ X \> X 2) linear - b _ a 



(88) 



[000164] The trouble arises when average in the orthogonal direction. This is what creates the 
troublesome antiderivaties. So in lieu of an exact average over an area, we use an approximation 
by computing the average of equation 88 evaluated at the endpoints of the domain & € [c, d]. 
This process is: 

1 



= W^aj^ b,C ^~ g ( a,C ) + 8 ( a ' d ^ 
[000165] Defining the vertex interaction matrix T\ as: 



-1 0 0 0 

0 10 0 

0 0 10 

0 0 0 -1 



and the lenslet width vector h as: 



h = 



b- a 
d- c 



allows us to write the approximation as: 



where the vertex vector V is given by: 

V T = ((a,c),(b,c),(b,d),(a,d)) . 



(89) 



(90) 



(91) 



(92) 



(93) 
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[0001 66] Similar machinations yield the result: 



(94) 



where the interaction matrix is: 



-10 0 0 
0 10 0 
0 0 10 
0 0 0 -1 



(95) 



[000167] The approximations in equations 92 and 94 are now the surrogate values for the exact 
form of the averaged slopes. The appendix explores the validity of these approximations and 
shows that it many cases it is exact. 

[000168] At last we can replace equation 87 and pose the new merit function as 



2 N ( f \ 2 

s M - — r v(v ) 



(96) 



liberating us from the problems of having to deal with antiderivatives. The minimization is 
routine. The quantities being varied are the wavefront values at each vertex. All the derivatives 
with respect to the wavefront values at the vertices are set equal to zero: 

*/ = 0. (97) 

[000169] The vector p now indexes the vertices, not the lenslets. 

[000170] The net effect is a series of linear equations that look like this (after the simplification 
hl=h2 = h): 
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«,v " Y „.,v-. ~ V „ + ,.v-. " V u+ ,v +I - V u -,.v + . = J (& M + ^(„ +1 , v) - s T 2 S {u v+l) - s(S (u+X v+l) ) (98) 
where s is an alignment spinor given by: 



s = 



1 
1 
1 

-1 



[000171] The vertex interaction matrix for this is: 



0 0 0 



(99) 



0 s, 0 0 



0 0 



0 



0 0 0 5, 



Ij 



(100) 



[000172] were 0 has the correct dimension, i.e. : 



0 = 



(101) 



[000173] Equation 98 is quite revealing in the way it relates the wavefront values y to the focal 
spot shifts 8. First, the vertex action is between next-to-nearest neighbors, not nearest neighbors. 
The nearest neighbor interactions cancelled out. Next, the focal spot shift content is quite rich, 
using the measurements of the four lenslets bound by the next-to-nearest neighbors. Finally, 
these equations are linear, so we expect a direct solution without iteration. 
[000174] This formula describes the special case for a lenslet on the interior of the array. These 
lenslet all have next-to-nearest neighbors. But many lenslets are on the edge and some are at the 
corners. The question is how should special shapes be handled? The answer is simple and 
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comes from considering the symmetry of the problem. When a photon hits a lenslet, the photon 
can't determine the pedigree of the lenslet. The photon is deflected and impacts the CCD array. 
The amount of the deflection is the same regardless of where the lenslet is located within the 
array. So we are guided to look for a solution which respects this symmetry. 
[000175] The answer is to include an existence function e u> v which tells us whether or not a 
lenslet exists. This also allows the user to construct arbitrary apertures. If the lenslet (u, v) 
exists, then e Ui v = 1, otherwise e Ut v = 0. Now equation 98 can be generalized to solve for 
arbitrary masks and apertures as sculpted by the existence function: 

^ ( r c T c T T \ 

= / P M S M + S2 >l/(«+l,v) " S 2% mV+ \) € (u,v+\) " 5 1 S (u + \,v + \) S (u + l, V+ \) ) ( 102 ) 

[000176] This is a series of + l)(r| + 1) linear equations which can be written in matrix form 
as: 

AV = D. (103) 
[000177] Now A is the vertex interaction matrix which describes the active vertices, v|/ is the 
vector containing the wavefront values at every vertex, and D is a vector containing 
combinations of the measured focal spot shifts. 

[000178] The matrix A is singular. A is singular if there is some nullspace x that maps A to 
zero. That is there exists some vector x such that A-x = 0. A look at equation 102 shows that 
each of the rows always sums to zero. Hence all constant vectors k are a nullspace for A. This 
means that we cannot find some inverse matrix A' 1 such that 
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A-'A^l. (104) 
[000179] But this does not mean that we cannot solve equation 103 with y. We are simply 
forced to compute a pseudoinverse A*. For linear least squares fit problems, the method of 
choice to form the pseudoinverse is the singular value decomposition, SVD. 
[000180] The process involves decomposing A into three matrices: a column orthogonal matrix 
U, an orthogonal matrix Vmd a diagonal matrix Z. The diagonal elements of 2 are the singular 
values of A. The process of SVD solves for these three matrices such that 

A=U-l-V T . (105) 
[000181] The pseudoinverse is now given by 

A + = V-1- X 'U T . (106) 
[000182] The entire process of reconstruction distills down to the matrix multiplication. 

V = A + Q . (107) 
[000183] This is a very powerful statement for it shows how most of the computation can be 
done in advance. For a fixed lenslet geometry, the matrix A* can be precomputed and stored. It 
does not need to be created during the measurement process. 

[000184] There is some fine print however. First of all, A + is very large. For a 64 x 64 lenslet 
array, has 65 4 elements which is about 136 MB of data. Next, the computation may take 
hours. So it behooves the user find ways to avoid recomputation when the aperture changes. 
Some strategies involve storing A + for many common apertures such as a series of circular 
masks. Another strategy is to pad the data. Consider A + for a full aperture. All lenslets are 
active. Then for some reason a cluster of three lenslets is excluded from the data. A good idea is 
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to use the same A* matrix and just pad the data in the excluded regions with some form of 

average data from the neighbors and exclude these vertices in the final result. 

[000185] So in closing, the zonal reconstruction starts with the antiderivatives and uses matrix 

multiplication to compute the wavefront values at each vertex in one step. 

[000186] In a very influential paper published in 1980, W. Southwell set forth what has become 

the standard method for wavefront reconstruction. The paper very clearly describes his 

reconstructor and is worth reading. Surprisingly though, Southwell solved his linear equations 

using successive overrelaxation (SOR). While SOR is an extremely valuable method, we feel 

that SVD is usually the better choice. 

[000187] We will compare the philosophical differences between the two reconstructors in the 
following section. Here we will compare the mathematics. If we define a neighbor spinor a as: 



which is just a collection of the basis vectors, then the Southwell merit function takes the form: 



where we are now dealing only with lenslets, not vertices. In this implementation, Southwell 
places the computed wavefront value in the center of the lenslet, not in the lower left- or upper 
right-hand corner as expected. 



a = 



0 

o : 



(108) 




(109) 
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[000188] The explanation for the placement of the wavefront value is based on the Southwell's 
treatment of the focal spot shift as relating to the wavefront slope - not the average slope. 
[000189] In the continuum, derivatives are measured at a point and are defined through the 
limit 



WW- r 



K ■ (110) 

[000190] In the discrete space of the lenslet h\ never reaches zero. It is limited by the size of 
the lenslet. In this space, the derivative is using the difference of neighboring points. Extending 
the discrete derivative is trivial. For the lenslet with a size vector given by equation 91, the 
partial derivative with respect to both variables is written 

V\,2g\ X \> X 2 ) * ^ • (HI) 

[000191] In light of this equation 68 can be viewed as 

wX x ^* d A x ^ x X 0*2) 

[000192] This seems to suggest that if the measurements are treated as slope measurements, 
and not averaged slope measurements, then the wavefront values should be placed in the lower 
left-hand corner of the lenslet or the upper right-hand corner depending on whether the limit is 
approaching zero from the right or the left. 

[000193] Returning to the Southwell minimization problem, we proceed as before. The 
parameters being computed are the wavefront values at the center of the lenslet and we vary A 
with respect to these values and set the derivative equal to zero: 
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^/ = 0. (113) 

[000194] For a £ x r| lenslet array with h\ = A 2 = A, this leads to N = ^r| linear equations of the 
form 

[000195] (The existence function has been eliminated to facilitate comparison). So the 
Southwell interaction involves nearest-neighbor lenslets and uses only four measured values. 
[000196] As before we can write a matrix expression just like equation 107. The lenslet 
interaction matrix A here is singular and we again resort to SVD. This allows for a direct 
solution and in many cases, we may precompute A + . Both methods will execute in equivalent 
times. 

[000197] It is time to step away from the mathematical details and talk about the differences 
between the two reconstruction strategies. The new method strictly adheres to the assumptions 
of Shack-Hartmann operation. It treats the measurements as being extended quantities, not point 
measurements. This is why the wavefront value cannot be placed at a single point. In this 
model, we are only permitted to make statements about the values at the vertices. Subsequent 
mathematics are then dictated. 

[000198] In Southwell's venerable model, he assumes that the difference between the 
wavefront values in adjacent lenslets is related to the average slope across those two lenselts. On 
thing that seemed unusual was the average of the two measurements was not placed on the 
boundary between the two lenslets. Regardless, the Southwell model looks at a lenslet, it 
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averages the 5i measurement with the lenslet in the 1 direction and averages the 5 2 measurement 
with the lenslet in the 2 direction. 

[000199] Southwell's model mixes information from lenslets immediately. The new model 
assumes that the measurement of a focal spot shift reveals information about the differences in 
wavefront values at the vertices. However, the neighboring lenslets each share two vertices. So 
the competing values for the shared vertices are adjudicated by the least squares fit. So clearly, 
the new model is also sharing information between lenslets, but not until the least squares fit. 
[000200] Perhaps the most obvious shortcoming with an immediate averaging is that small 
features of the order of a lenslet will be diminished. In day to day practice, this will probably not 
be noticed by most users. Only those requiring exquisite detail resolution would be driven 
towards the new method. 

[000201] Another interesting difference between the two is in the wavefront value interaction. 
The Southwell equations balance the wavefront values with the four nearest neighbors. In the 
new method, the nearest neighbors interactions cancel out, and the surviving interaction between 
the next-to-nearest neighbor vertices. This is shown schematically in FIGs. 14-15. FIGs. 14-15 
compare the interactions of the wavefront values in the least squares fit. In the Southwell model 
depicted in FIG. 14, the wavefront value is balanced between a central lenslet (marked in red) 
and the four nearest neighbors. This can be seen in equation 1 14. In the new model shown in 
FIG. 15, the wavefront value is balanced between a central vertex and the four next-to-nearest 
neighbors and the nearest neighbor interactions exactly cancel. The resultant interaction is 
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specified in equation 98. So it is not that the nearest neighbor interaction is avoided, it is just that 
they exactly cancel out. 

[000202] Finally we note that the measurement content of the 0 matrix is richer in the new 
model. Southwell blends together four measurements, and the new model blends eight. 
[000203] In summary, both methods have a comparable run time and the memory requirements 
are virtually identical. They are both direct methods and both allow for arbitrary aperture shapes. 
The primary difference seems to be an improvement in resolution of fine details. Comparisons 
in the way the two models handle noise need to be done yet. 

[000204] This section evaluates some sample wavefront that we put into both models and give 
the reader a visual comparison of the two strategies. The input wavefronts are shown in FIGs. 
16-17. FIGs. 16-17 shows two different wavefronts to highlight the differences between the two 
reconstructors. The first case in FIG. 16 presents three bumps with different sizes. This shows 
how the two methods are able to resolve features of different sizes. The second case in FIG. 1 7 
was created to be difficult or impossible to reconstruct. The results were surprising. 
[000205] FIG. 16 shows a collection of three Gaussian bumps. In practice, reconstructors have 
trouble picking out small features against a background of big features. This set of bumps was 
designed to probe this problem. The second wavefront resembles an egg crate and it was created 
to defy reconstruction and examine the behavior of the reconstructor as it fails. 
[000206] Turning to the first case in FIG. 16, the goal was to challenge the reconstructor in an 
attempt to resolve smaller bumps in the presence of a large bump. Here we see two rather sharp 
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bumps in front of a prominent feature. The smaller bump on the left is taller and wider than the 
bump to the immediate right. This bump is the most difficult one to resolve. 
[000207] The results are shown in FIGs. 18A-19 and FIGs 20-21B, comparing the new 
reconstructor to the Southwell reconstructor. The wavefronts shown in FIGs. 16-17 were 
submitted to the two reconstructors. FIG. 19 shows the input wavefront at the resolution of the 
lenslet array and FIGs. 18A-B show the results of the two reconstructors. The new method 
recreates the smaller bumps better indicating that this method is best suited to resolving fine 
details. FIGs 20-2 IB show a series of two dimensional plots because at the resolution of the 
lenslets, the three dimensional plots are difficult to interpret. Suprisingly, the output from the 
new reconstructor (FIG. 21 A) matches very well with the input. The Southwell method (Fig. 
2 IB) did not fare as well, but still performs adequately. This suggests that the new reconstructor 
is better suited to handled complicated structure. 

[000208] Since the smaller bumps are the features of interest, the plots have been scaled to their 
height. The input wavefront of FIG. 19 has been plotted at the same resolution as the lenslet 
array to facilitate comparison with the reconstructor outputs. As expected, the new reconstructor 
of FIG. 18B fared better at resolving the small bumps. The bump on the left is faithfully 
reproduced. The bump on the right is reproduced vaguely. All reconstructors have an amplitude 
ambiguity, so we quantified the problem by essentially doing at least squares fit to see what 
magnification it would take to get the reconstruction to most closely match the input. The new 
reconstructor overstated the amplitude by 0.4%; the Southwell reconstructor understates the 



55 



Docket No. WFS.017 
-- PATENT - 

amplitude by 2.7%. The x 2 values were 0.003 for the new method and 1 1 .6 for the Southwell 
method. 

[000209] After running dozens of cases, we wanted to present a wavefront that could not be 
reconstructed in order to study the behavior as the algorithm failed. We wanted to identify 
artifacts associated with an errant reconstruction. We created the wavefront shown in FIG. 17. 
The magnification computation shows that the new reconstructor overstates the amplitude by 
11% and that the Southwell method understates the amplitude by 19%. The % 2 values were 
shocking. The new method has a x 2 of 2 10" 25 ; this represents the limit of double precision 
rational. In other words, the new method had a x 2 of zero. For the Southwell method the value is 
11,800. 

[000210] The fact that the new reconstructor had essentially no fitting errors was astounding 
and led to a careful exploration of the approximations in equations 92 and 94. The results, 
detailed in the appendix, show that the approximation allows for a very rich wavefront structure 
across each lenslet. It probably allows for more structure that the focal spot location method can 
tolerate. That is, the focal spot will degrade and become ambiguous before the approximation 
used in the reconstructor breaks down. 

[000211] This paper details new reconstruction strategies formally based upon the principles of 
the Shack-Hartmann sensor. This method was compared to a modernized implementation of the 
standard method developed by Southwell. The difference between the two different procedures 
is negligible or extreme depending upon the application. For common measurements of 
smoothly varying wavefronts, the difference between the two methods teeters on the edge of 
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imperceptible. However, if the wavefront of interest is richly structured or finely detailed, then 
the new method is a substantial improvement. 

[000212] The results on the egg crate reconstruction FIG. 17 were startling. After all, the 
wavefront was designed to break the reconstructor, yet the % 2 was zero to the limit of machine 
precision. This inspired a detailed look at the approximation used in equations 92 and 94. This 
section documents the investigation. 

[000213] To recap, the merit function (equation 87) for Shack-Hartmann sensors involves 
antiderivatives, which are unwelcome in a zonal reconstruction. An approximation was used to 
bypass the nettlesome antiderivatives and formulate a merit function based on the wavefront 
values. The approximation is contained in the two equations: 

[000214] We start by forming the sum and the difference of the two equations represented 
above. Remembering that the domain is x\ e [a, b], x 2 e [c, d\ we find that: 

T(M)-?M^DT 2 y (116) 

and 

^(b^-^ia^^D^^x^jD^^x,). (117) 

[000215] At this point we need to be more specific about y(jti, x 2 ). We can represent this 
function as a Taylor series. This is not because we champion the use of the Taylor monomials in 
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modal reconstruction; it is because any piecewise continuous function can be described exactly 
by an infinite series of monomials. We now write 

/=0 j=0 

[000216] When this expansion is substituted into equations 1 16 and 1 17, we see: 
ii^-jM J d J -a^) 

i=0 y=0 

-Hi 777 (* w - « w x^' -^'Vjii rft-Ab'-'" - a-»<id> - c'), (H9) 

«2 /=0 y=0 7 + 1 "i ,=0 7=0 * ~ 7 + 1 



and 



i=0 7=0 



■ rZ Z t-t(*'" j - « W X" / " - fZ Z t^(» W4 ' - "-"X* - c). (120) 

"2 /=0 7=0 7 + 1 "l i=0 7=0 2 " J 

[00021 7] Using the identity: 

■=2y~v ( 121 ) 



1=0 

one can show that the approximations in equation 1 15 are exact whenever the wavefront over the 
lenslet has a structure like: 

00 00 

V(x l ,x 2 )=a 00 + a u x x x 2 + £ a i0 x[ + £ a /0 *J . (122) 

[000218] This allows for an incredible amount of structure over the lenslet. For example, any 
separable function: 

58 



Docket No. WFS.017 
PATENT 



^(x l9 x 2 )^x(x x )Y(^) 



(123) 



will be handled exactly; in other words, there is an exact solution and it is described by equation 
107. Or any separable function with the lowest order torsion term, for example: 



which is exactly the form of equation 123. 

[000220] EFFICIENT COMPUTATION WIH SPECIAL FUNCTIONS LIKE THE CIRCLE 
POLYNOMIALS OF ZERNIKE. 

[000221] The circle polynomials of Zernike play a prominent role in optical analysis. While 
decompositions of wavefronts into Zernike polynomial series can yield valuable insight, 
computing with the polynomials themselves is quite inefficient. Here we outline how rational 
polynomials like those of Zernike, Legendre, Chebyshev and Laguerre can be handled as affine 
combinations of a Taylor monomial set. We demonstrate how calculations can be performed 
much more rapidly in the Taylor basis and how to use integer transformations to recover the 
exact amplitudes in the desired basis. We also explore C++ optimizations for storing the Zernike 
amplitudes and transforming between Zernike polynomials and Taylor monomials. 




(124) 



[000219] The egg crate wavefront in FIG. 1 7 is: 




(125) 
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[000222] The circle polynomials of Nobel laureate Friz Zernike have proven to be an 
invaluable tool for investigating optical systems. Zernike's polynomials are orthogonal which 
facilitates associating individual terms with physical observables. Generations of optical 
scientists have used this representation to describe optical system performance, to quantify 
behavior of individual components and to provide a framework for consistent evaluation of 
optical phenomena. All the quantification is done in terms of the amplitudes for each term in a 
Zernike polynomial expansion. While Zernike's polynomials have proven to be an invaluable 
mechanism for describing the physical world. The formulation is quite cumbersome and actual 
computations in this basis set are needlessly complicated. 

[000223] Because of the Gram-Schmidt orthogonalization process that Zernike used to 

construct his polynomials, the complexity of higher order terms grows very rapidly. Fortunately, 

the Taylor monomials form a minimal complexity representation that is exactly equivalent to the 

Zernike basis on an order-by-order basis. This paper formally establishes the equivalence of the 

two sets and then discusses how to compute in the Taylor basis and transform the results into the 

useful and familiar Zernike basis. In fact, this procedure is even more general: the Taylor 

monomials can be used to represent any rational polynomial set, such as those of Legendre, 

Laguerre and Chebyshev. With this realization, the Taylor basis becomes even more valuable 

for it is in essence all the other rational polynomial basis sets. So by computing in the Taylor 

basis, one is able to use affine transformations to recover the amplitudes of all other rational 

polynomial basis. This means, for example, that there needs to be only one reconstruction. With 

the Taylor amplitudes in hand, we are only one matrix multiplication away from the amplitudes 

of Zernike, Legendre, Chebyshev or Laguerre. And given the amplitudes in these other sets, a 
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matrix multiplication converts them into Taylor amplitudes, and the another matrix 
multiplication can take us into the basis set of choice. 

[000224] This application discuss some practical examples in wavefront sensing where the 
Zernike polynomials are cumbersome and tedious to use. We explain how to use the Taylor 
basis for computation and then how to recover the Zernike amplitudes. Besides laying out the 
fundamental mathematics, we go into to detail to show how this would be implemented in a 
modern computer language such as C++. Forethought in designing the architecture of the code 
results in an extremely facile and powerful tool for wavefront analysis. 

[000225] The mathematics represent some level of abstraction, and the computer code shows a 
definite interpretation. Several issues in implementation, such as run time and memory usage, 
are explored. Finally testing and debugging a large volume of code requires stringent validation. 
We outline the five layers of validation we used to verify several million lines of computer code. 
[000226] EQUIVALENCE OF THE ZERNIKE AND TAYLOR AMPLITUDES 
[000227] When the complete set of Zernike 6 and Taylor functions through arbitrary order d are 
considered, the two sets are equivalent. Consider the explicit example of two sets for second 
order. 



Taylor Function 



Zernike function 



T 2o( x >y) 
T n( x >y) 



xy 



y 



X 



1 




y 



x 



2xy 

2x 2 + 2y 2 -\ 

2 2 

y - x 



1 
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[000228] The equivalence is shown below. First, we express the Zernike functions in terms of 
Taylor functions: 

z oo( x >y)= T oo(x,y) 
z wi x >y)= T wi x >y) 

Z u (x,y) = T 0l (x,y) 
zJx t y)*2T n (x,y) 

Z 2l (x,y) = 2T 20 (x,y)+ 2T 02 (x,y)- T m (x,y) 
z 22( x >y) = T<n( x >y)- T m (x,y) 

[000229] Next we express the Taylor functions in terms of the Zernike functions: 

T oo( x *y) = Z oo( x >y) 
T l0 (x,y)= Z 10 (x,y) 

T m{*>y) s z n( x >y) 

T 20 {x,y)=^ z oo(x,y)+ jZ 2l (x,y)+ ^Z m (x,y) 

%n(*>y) s jZn{x>y)+ \ z n(x,y)+ \ z oo( x >y) 

[000230] So for example the arbitrary wavefront y(x,y), written as: 

* ( x >y) = too + 'io* + t m y + ' 2 o* 2 + tn*y + toiy 2 ( 12 6) 

can be expressed in terms of the Zernike polynomials as 



V( x >y)= ('00 + y + t -f)z m {x,y)+ t lo Z l0 (x,y)+ t m Z n (x,y) 



+ Y 2 ° x,y 



ft t 

l 02 l 02 

I 4 4 J 



Z 2l x,y + 



tp2 

V 2 



t 



\ 

20 
2J 



Z 22 x,y 



(127) 



62 



Docket No. WFS.017 
-PATENT- 



[000231] One already sees an inkling of the native simplicity of the Taylor basis. The 
relationship between the two bases can be cast in matrix form as 







Z l0 (x,y) 




z n {x,y) 




z 20 (x,y) 




z 2l (x,y) 




Z 22 (x,y) 





1 

0 



0 0 

1 0 
0 0 1 
0 0 0 
-10 0 



0 0 0 

0 0 0 

0 0 

2 0 

0 2 



0 0 0-101 



T oo{x>y) 

T \J<X,y) 

T oi(x>y) 
T 2o( x >y) 
T u(x>y) 
T <a{x>y). 



(128) 





l 


0 


0 


0 


0 


0 " 






Too(x,y) 




0 


1 


0 


0 


0 


0 


Z(X) 


[ x >y) 


T io(x,y) 




0 


0 


1 


0 


0 


0 


z l0 


[**y) 


T 0i {x,y) 




1 


0 


0 


0 


1 


1 


zj 


[x,y) 


T 20 (x,y) 




4 


1 


4 


~ 2 


z 20 


(x,y) 


T u {x,y) 




0 


0 


0 


2 


0 


0 


z 2i 




jjx,y) 




1 


0 


0 


0 


1 


1 


z 22 






.4 


4 


2 . 







(129) 



[000232] But even this notation is cumbersome and it is more convenient to write these 
relationships in the forms 

z(x,y)=fz T T{x,y), (130) 

or 

T(x,y)=tz T -Z{x,y). (131) 

[000233] The matrix tz transforms an amplitude vector from the Taylor basis to the Zernike 
basis and the matrix fz transforms an amplitude vector from the Zernike basis into the Taylor 
basis. Their usage in equations 128 and 129 may seem counter-intuitive, so a short 
demonstration follows. 
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[000234] The mappings in equations 128 and 129 show that the Zernike and Taylor basis are 
equivalent. This means that any wavefront \y(x, y) can be resolved in either basis allowing us to 
write either: 

v(x,y)=t T T(x 9 y) (132) 

or 

v(x,y)=z T -Z(x,y) (133) 

where t and z are respectively, amplitudes in the Taylor and Zernike spaces. These amplitudes 
are the outputs from modal reconstruction. Notice that the results of the dot products are a scalar 
function, in fact, they are exactly the same scalar function as seen in equation 127. Since 
equations 132 and 133 are the same function, we can write: 

z T .z(x,y)=t T -T(x,y) (134) 

which in light of equation 130 can be recast as: 

z T .(fz T .T(x,y))=t T -T(x,y). (135) 

[000235] This immediately implies: 

z T -fz T =t T (136) 

which is of course: 

t=fz-z. (137) 
[000236] A similar process yields the relationship: 

z=tz-t. (138) 
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[000237] The naming scheme for the transformation matrices should be now be clear. The 
matrix tz takes amplitudes from the Taylor basis to the Zernike basis and fz takes amplitudes 
from the Zernike basis into the Taylor basis. As one might expect, these operations are the 
inverse of one and other. That is: 

tzfz = fz-tz= 1 (139) 
where 1 is the identity matrix of appropriate order. 

[000238] This analysis can be extended arbitrarily order by order. The approximation theory of 
Weierstrass, p. 1929 shows that in the limit d either set of polynomials can be exactly 
represented any piece wise continuous wavefront. 
[000239] The matrix fz is based upon the recursion relationship: 



q m m- J 



Z nm(x>y)= X Z Z (~ l T I 9 - n I 7 I —( XT, AT 



2{i+k) +Py n-2(i+j+k)-p 



(140) 



where s= sgn(v),p= ^{s+ l),q= ^-(v- sMod 2 n), and M = 



n 

m-- 



. For a given order n and 



index w, this formula can be used to generate a vector of coefficients for the fz matrix. When the 
fz matrix is fully assembled, the tz matrix is computed as 

tz=fz\ (141) 
[000240] Of course the transformations represented,^ and tz are not just a single matrix; they 
are a series of matrices, one pair for each order. 

[000241] By comparison, the Taylor monomials are generated by a much simpler function: 

T i _ JJ (x,y)=x i - J y J . , (142) 
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[000242] In summary, it is perhaps convenient to consider the two bases, Taylor and Zernike, 
as having a relationship similar to the one between the Cartesian and spherical polar coordinates. 
Certainly we all realize that (x t y, z) references the exact same place as (r, 9, OJ. These 
mappings are both complete and unique. They are completely equivalent ways to express 
location. Integrations in both coordinates yield identical results. 

[000243] We are all familiar with the Taylor expansion and we know that any piecewise 
continuous function can be exactly reproduced by an infinite number of expansion terms. The 
Taylor monomials represent the simplest basis for an expansion. All other expansion sets, such 
as those of Laguerre or Legendre, are just combinations of the Taylor monomials. The natural 
question then is why are these other bases so popular? It is primarily due to the fact that the 
expansion polynomials are orthogonal. In general, two functions fi(x,y) and fi(x,y) are orthogonal 
over some domain D with respect to a weighting function w(x,y) if they satisfy: 

J lfi{x 9 y)fi(x 9 y)^x 9 y)dxdy = (143) 

D 

where oty is some constant and 8,y is the Kronecker delta tensor. 

[000244] We are all familiar with the concept of orthogonality since we live in a Euclidean 
world. Consider some point P on an x-y graph. We are able to measure x and y independently 
because the axes are orthogonal. As we vary in the x- direction at P, the y values do not change. 
In general for a coordinate system we can formulate the geodesic between two points by 
integrating ds using the metric g: 

ds 2 = YjSijdqA 2 (144) 

u 

[000245] For the Euclidean world that we live in, dq\ = dx, dqi = dy, and the metric is: 
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1 0 

'-[o ij- ' (145) 

[000246] In this special case, dx and dy may have finite size. In general, orthogonal coordinate 
systems have gy = 0 when i j: 

gij = tyy (H6) 

where the by are constants or functions. This states mathematically what our intuition tells us: 
we can be at a point P and vary our position in one coordinate without affecting the values in the 
other coordinates. 

[000247] A comparison of equations 143 and 145 suggests the concept of orthogonality extends 
in a natural way to functions. The consequence is that measurements of the amplitudes are 
independent of each other. For example, consider the Zernike polynomials. It is possible to do a 
least squares fit to measure a single amplitude like Z42. The measurements are independent and a 
computation of subsequent amplitudes will not affect the previously measured amplitudes. On 
the other hand, a non-orthogonal set of basis function like the Taylor monomials behaves 
differently. The amplitude form a term like x 3 y 2 depends upon which terms are being measured 
with it. This makes the matter of physical interpretation quite difficult. Consider the lowest 
order Zernike polynomial term, Zoo. It's amplitude is the mean value of the surface being 
analyzed. The corresponding Taylor amplitude - the constant term - will fluctuate depending 
upon how many terms are in the fit. 

[000248] But as we will see, it is easier to do computations in the Taylor bases and then 

transform the amplitudes into an orthogonal set. The Taylor monomials are the simplets bases 

set to work with. FIG. 23 shows schematically how the functional bases are related. Think of 
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the Taylor monomials as a hub connecting all the rational polynomials basis sets; four such sets 
are shown here. With software designed around the Taylor basis, the user assures the fastest 
possible execution and the least amount of computer code. An affine toolkit is used to transform 
the results of the Taylor-based calculations into the other basis sets. These affine 
transformations are integer based and virtually instantaneous. The Taylor basis is the hub and is 
the set of minimal complexity. After doing computations in the Taylor basis, rational 
transformations are used to compute quantities any desired basis. This also allows for rapid 
conversion between any two bases such as Zernike as Chebyshev (with the usual caveat that the 
domains must match). For example, if/c and cz are the transformations that take amplitudes 
from and to the Chebyshev basis, then we could start with a vector of Zernike amplitudes, z, and 
compute the vector of Chebyshev amplitudes, c, via: 

c =tc-fz-z. (147) 
[000249] For the case of n polynomial bases, there are 2n transformations. If we don't use the 
Taylor bases as the hub and have transformations directly connecting all the bases, then n(nA)/2 
transformations are needed. The Taylor hub will greatly simplify the affine transformations. 
[000250] FIGs. 24A-B illustrates that aperture resizing is just rescaling the coordinate system. 
Consider an arbitrary function shown in FIG. 24A). We can shrink the x-coordinate by a factor 
of/= Vi, as shown in FIG. 24 B. The curve does not change. We have effectively just relabeled 
the coordinate system. 

*W= J™ZS«/-yX 7 - ( 147 ) 
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[000251] Of course practical applications truncate the expansion at some finite d. Now we 
apply the transformation in equation 146 to 147 and we see that the amplitudes of order d are 
multiplied by the factor/. For example the transformation matrix F that adjusts the amplitude 
vector for a second order expansion is 



F = 



f° 


0 


0 


0 


0 


0 


0 


/' 


0 


0 


0 


0 


0 


0 


/* 


0 


0 


0 


0 


0 


0 


f 2 


0 


0 


0 


0 


0 


0 


f 2 


0 


0 


0 


0 


0 


0 


f 2 



(148) 



[000252] Now if the user had Zernike amplitudes foisted upon him and wished to perform the 
resize operation, the affine toolkit will help. Given a Zernike amplitude vector z and a resize 
matrix F, the new Zernike amplitudes z would be given by 

z' = tz-F-fz-z . (149) 
[000253] And were we given Zernike amplitudes and asked to reduce the aperture and calculate 
the Chebyshev amplitude vector c we would compute 

c-tcFfzz, (150) 
[000254] Yes, it may be possible to compute F in other bases, but it is trivial in the Taylor 
basis. And working in this basis allows the user to compute amplitudes in any other rational 
polynomial basis set. 

[000255] While in general the Taylor polynomials are the most efficient computation basis, the 
rule is not absolute. There are instances when the evaluation is much faster in the Zernike basis. 
These cases involve wavefronts with rotational symmetry, such as spherical waves. 
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[000256] To study how increasing the fit order d improved the resolution of sharp details, a 
cone was modeled. The equation for a cone in polar coordinates is: 

*(r,0) = 1-r. (151) 
[000257] The strategy is to perform a least squares fit to find the amplitudes which best 
describe the cone. If we compute the Zernike amplitudes, then we may compute them one order 
at a time. The orthogonality of the Zernike polynomials guarantees that as higher order 
amplitudes are computed, the lower order amplitudes will not change. This is the great 
advantage of an orthogonal basis set. The merit function for order n is written as 



Z 2 = \\[^{r,e)-tc nm Z nm (rA rdrdO. 



(152) 



[000258] The merit function is minimized by setting derivatives with respect to the amplitudes 
equal to zero and solving the resultant set of linear equations. For example, picking the arbitrary 
amplitude eg: 



V 2 = 0 - 

[000259] This immediately leads to the equation: 

l][t c m ZjrAz u (r,e)rdrd6= \\v(r,9)Z,{r,8)rdrde 



(153) 



o 0 V m=Q 



In 1 



(154) 



0 0 



[000260] There is very little integration needed to solve this problem. The left hand side is 
solved by using the orthogonality relationship: 



2- 



sgn 



n 

m- — 



2J) 



(155) 
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and we see that this integral projects out the value of Cy. The right hand side is greatly simplified 
because of the periodic nature of trigonometric functions. Since sin (0 + 2ri) = sin (0), the only 
integrals that can be non-zero are the rotationally invariant terms which do not have a 0 
dependence. The general solution for the only surviving coefficient is given by 



n\ 1 



\\v(rfydrdO. (156) 



"•2 * 0 0 

/ 

[000261] This is greatly simplifies the problem. In a few minutes we were able to perform a 

Zernike polynomial fit to order 1000. This expansion series looks like: 

[000262] 

1 _ 2 2 2 _2 2 2 2 2 2 

* = 3 Z 0,0 " 5 Z 2,l + 21 Z 4,2 " 45 Z 6,3 + 77 Z 8,4 " j j ? Z ,0,5 + ^ Z 12,6 " ^ Z I4,7 + ^ Z 16,8 " ^ Z 18,9 

2 2 2 2 



437 Z 2o,io 525 Z 2o,io + - + 10 i97 Zl00 ' 50 *" + 1001997 Zl000 ^ 500 * (15?) 
[000263] Attempting this feat in the Taylor basis would require the inversion of a 501,501 x 
501,501 matrix, so this method was much faster. The issue here is that this problem is analytic, 
not numeric. Applications involving measurements are typically numerical, and in those cases 
the Taylor basis is much faster and easier. 

[000264] FIG. 25 shows some results for this fit, using the Zernike basis to speed computation. 
The first image shows the input function in equation 151. Subsequent images show the fits to 
successively higher orders. The first order term, the piston, is the mean value of the function. 
The higher order terms add curvature. Eventually they will sharpen the point and straighten the 
sides. By order 40, the fit closely resembles the input curve. Starting with a rotationally 
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invariant object, the cone shown in the upper left-hand corner, successive least squares fits were 
done to extract the amplitudes of the rotationally invariant Zernike polynomials. A fit to order 
1000 took only a few minutes. This computation may not be feasible using the Taylor basis. 
The results of fits to various order are shown in subsequent frames. 

* 

[000265] At this point, we hope to have established the background for and strongly motivated 
the use of the Taylor monomial basis set. The final component is to discuss the issues of putting 
these ideas into an efficient, rapid and robust computer code. Some forethought is required, but 
after starting with a good design, the code is very easy to use. 

[000266] To support the aforementioned methodologies of computation, data structures and 
algorithms need to be optimized for this solution. Because of the flexibility and speed 
requirements of these computations this implementation is done in C++. 
[000267] A three dimensional array can be used for the representation of the Zernike 
polynomials. By using a three dimensional array, it becomes simple to access any defined order, 
polynomial or coefficient. This first dimension of the array is the order and is addressed Z[o]. 
The second dimension of the array is the polynomial and is addressed Z[o][p]. Finally the third 
dimension of the array is a coefficient, which is addressed Z[o][p][c]. 

[000268] Using a three dimensional array can be very inefficient if not implemented correctly. 

As orders increase, so do the number of polynomials and coefficients. To maintain our 

addressing structure, a statically allocated array would be very wasteful. Consider building the 

Zernike array to third order. A statically allocated array would have to be declared INT64 

Z[3][3][10]. This array statically allocates 3 polynomials and 10 coefficients for all orders even 

though orders 0,1, and 2 have fewer polynomials and coefficients. By using dynamically 
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allocated memory we can prevent ourselves from using excess memory and still maintain our 
addressing structure. Building the Zernike structure up to but not including n-Order is actually 
very trivial. Listing 1 below shows a sample allocation. 
[000269] 



typedef INT64 Coeff; 

typedef Coeff* Poly; 

typedef Poly* Order; 

typedef Order* zernikes; 



// Coeffs 64-bit integers 

// Polys array of Coeffs 

// Orders array of Polys 

// zernikes array of Orders 



int n Poly; 
int nCoeffs; 



// Num poly 
// Num coeffs 



try{ 

z = new Order [nOrder] ; 
for( int i-0; i<nOrder; i++ ) { 
nPoly = i + 1; 

nCoeffs = ( (i+1) * (i+2) ) /2; 

z[i] = new Poly[nPoly]; 

for (int j=0; j<nPoly; j++ ) { 

z[i] [j] = new Coef f [nCoeffs] ; 

memset (z [i] [j] , 0, size of (Coef f ) * (i+1) ) ; 



} 



} 



} catch ( std: -.bad alloc ) { 
FreeZernikes (z) 
Return; 

} 



// Try to allocate 

// Allocate orders 

// For each order 

// Polynomial in order 

// Coeffs in polynomial 

// Allocate polynomials 

// For each polynomial 

// Allocate coeffs 

// Set coeffs to zero 

// Matches inner for 

// Matches outer for 

// Catch alloc error 

// Free memory 

// Exit function 

// Matches catch 
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[000270] Constructing the Zernike coefficient table in this fashion avoids unnecessary 
allocation of memory. After constructing this structure, the non-zero coefficients need to be set 
to the correct value. This can be done in one of two ways. First by calculating the non-zero 
coefficients for every polynomial in all orders. The following code segment generates these non- 
zero coefficients for any Zernike polynomial where n is the order and m is the offset of the 
Zernike in order n. 
[000271] 

Int64 c, p, d, S; 
Double M, q; 

long mu, nu, kappa; // index variables 



d = n - 2*m; 
s = Sgn(d) ; 

M = (n/2.0) - fabs(m -(n/2.0)); 
P = (abs(s)/2.0) * (s + 1) ; 
q = s * (d - s*(n%2) ) / 2.0; 



// Angular frequency 

// Side of periodic table 

// Dist from center line 

// cos/sin term determinant 

// Loop variable 



for (long i = 0; I <= q; i++) { 
for (long j = 0; j <= M; { 

for(long k = 0; k <= M- j ; k++) { 
mu = 2 * (i+k) + p; 
nu = n = 2 * (i+j+k) - p; 
kappa = mu + nu; 

c = pow(-l # i+j) * Binomial (n-2*M, (2*l)+p) * 
Binomial (M-j,k) * (nFact (n-j ) / (nFact ( j ) * 
nFact (M-j) * nFact (n-M-j) )) ; 

INT64 index = (kappa * (kappa+l))/2 + nu; 
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// Evaluate a coeff 
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Z[n] [m] [index] += c; 



} 



// Update the coeff 

// Matches inner for 

// Matches middle for 

// Matches outer for 



[000272] It's obvious that this calculation is at least a ®(n 3 ) operation due to the multiple 
embedded for loops. The second method is to load the non-zero amplitudes into the array 
manually, the previous section of code can be easily modified to generate the C++ code to make 
this a 0(1) operation. 
[000273] 







II 


0 th order 


coeff 


z[0] [0] [0] = 


(INT64) 1; 


II 


n = 0, m 


= 0 ; constant 






II 


1 st order 


coef f s 


z[l] [0] [1] = 


(INT64) 1; 


II 


n = 1, m 


= 0 ; x 


Z[l] [1] [2] = 


(INT64) 1; 


II 


n = 1, m 


= 1; y 






II 


2 nd order 


coef f s 


z[2] [0] [4] = 


(INT64) 2; 


II 


n = 2 , m 


= 0; xy 


z[2] [1] [0] = 


- (INT64) 1; 


II 


n = 2, m 


= 1; constant 


z[2] [1] [3] = 


(INT64) 2; 


II 


n = 2, m 


= 1; x 2 


z[2] [1] [5] = 


(INT64)2; 


II 


n = 2 , m 


= l; y 2 



z[2] [2] [3] = -(INT64)1; 
z [2] [2] [5] = (INT64) 1/ 



// n = 2, m = 2; x 2 
// n = 2, m = 2; y 2 



75 



Docket No. WFS.017 
—PATENT— 



// 3 rd order coeffs 



z [3] [0] [6] 




- (INT64) 1 ; 


II 


n 




3, 


m 




0; 


x 3 


z[3] [0] [8] 


- 


(INT64) 3; 


II 


n 


= 


3, 


m 


_ 


0; 


xy 2 


z[3] [1] [1] 


= 


- (INT64)2; 


II 


n 


= 


3, 


m 


= 


1; 


X 


z [3] [1] [6] 




(INT64) 3 ; 


II 


n 




3, 


m 




l; 


x 3 


z[3] [1] [8] 


- 


(INT64) 3/ 


II 


n 


_ 


3, 


m 


_ 


1; 


xy 2 


z[3] [2] [2] 




- (INT64) 2; 


II 


n 




3, 


m 




2; 


y 


z[3] [2] [7] 




(INT64)3; 


II 


n 




3, 


m 




2; 


x 2 y 


z[3] [2] [9] 




(INT64) 3; 


II 


n 




3, 


m 




2; 


y 3 


z[3] [3] [7] 




- (INT64) 3; 


II 


n 




3, 


m 




3; 


x 2 y 


z[3] [3] [9] 




(INT64) 1; 


II 


n 




3, 


m 




3; 


y 3 



// 4 th - nth order coeffs 

[000274] Since this methodology uses only assignment operations, and lacks any mathematical 

computations or for loops, this segment executes in a constant time. The second method makes a 

larger data segment of an executable, DLL or a statically linked library. To illustrate the 

differences, a sample program that builds and frees the Zernike array 1000 times was created and 

the coefficients were populated using both methods. The results are summarized in table 1. The 

first method creates an executable that is 168KB. The same sample program using the second 

method creates an executable that is 2 MB. When these two programs were ran on a machine 

with two 1.5 MHZ processors, the first method took on average 385 ms while the second method 
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only took 14 ms on average. Clearly speed is achieved by using more disk space. When speed is 
an issue the second method should be used. 

[000275] One other issue comes into play while deciding which method to implement. The 
coefficients that represent the Zernike polynomials through order 39 are small enough to be 
represented by a 64-bit integer. However, due to the use of factorials in the computation of these 
values, an arbitrary precision number type is necessary to compute them. In fact, using a naive 
implementation without arbitrary precision, the maximum order of the coefficients that can be 
correctly computed is 20. To calculate the coefficients during the runtime, ref.[3] outlines 
functions necessary to write a C implementation that handles extremely large numbers and the 
algorithms they chose are highly efficient. Otherwise, a tool such as Mathematica will prove 
itself priceless for generation of the coefficients, since all of its calculations are done using 
arbitrarily precise numbers. It is also necessary to note that above order 40, the coefficients 
become large enough to overflow a 64-bit integer data type. 
[000276] Table 2: Methods for Loading Amplitudes 



Method 


Time to Load 
Amplitudes 


Memory Used 
During 
Execution 


Memory Used 
on Disk 


Computing 


385 ms 


1388 KB 


168KB 


Amplitudes 








Pre 


14ms 


3312 KB 


2MB 


Computed 








Amplitudes 









Now that we have discussed the internal representation of the Zernikes, the computation 

transforming them into their Taylor representation needs to be considered. The matrix 

operations in equations 128 and 129 are an example of how the transformations are performed 
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for order 2. Matrix multiplication quickly becomes very expensive. For example the 
transformation matrices at order 30 are 496 x 496. However, since the matrices are very sparse, 
it is extremely time and memory inefficient to create the transformation matrices and then 
perform a multiplication. Using a similar technique as the one used to create the Zernike vectors, 
we use equation 140 to create code to perform sparse matrix multiplication. For example, 
converting 4 th order Zernike polynomial to its Taylor equivalent can be performed as follows: 
[000277] 

tay[0] = zer[0] + zer[3] / (INT64)4 + zer[5] / (INT64)4; 
tay[l] = zer[l] + zer[6] / (INT64)2 + zer[8] / (INT64)6; 
tay[2] = zer[2] + zer[7] / (INT64)2 + zer[8] / (INT64)2; 
tay[3] = zer[4] / (INT64)2; 

tay[4] = zer[3] / (INT64)4 + zer[5] / (INT64)4; 
tay[5] = -zer[3] / (INT64)2 + zer[5] / (INT64)2; 
tay[6] = -zer[6] / (INT64)4 + zer[8] / (INT64)4; 
tay[7] = zer[6] / (INT64)4 + zer[8] / (INT64)12; 
tay[8] = zer[7] / (INT64)12 + zer [9] / (INT64)4; 
tay[9] = -zer [7] / (INT64)4 + zer [9] / (INT64)4; 

where tay is the array holding the Taylor coefficients and zer is the array of Zernike 
coefficients. With this implementation, the transformation of the Zernike coefficients to the 
Taylor representations performs only 19 divisions, whereas the matrix form would do 225 
divisions. To do the inverse all that is necessary is to perform these operations: 
[000278] 

zertO] = tay[0] - tay [4]; 

zer[l] = tay[l] - tay [7] * (INT64)2; 
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zer [2] 




tay[2] 




tay [8] * 


(INT64)2; 




zer [3] 




tay [4] 


* 


(INT64) 2 


- tay [5] ; 




zer [4] 


= 


tay[3] 


* 


(INT64)2; 






zer [5] 




tay [4] 


★ 


(INT64)2 


+ tay [5] ; 




zer [6] 




-tay [6] 




* tay [7] * 


(INT64)3; 




zer [7] 




tay [8] 


★ 


(INT64) 3 


- tay [9] * 


(INT64) 3; 


zer [8] 




tay [6] 


* 


(INT64) 3 


+ tay [7] * 


(INT64) 3; 


zer [9] 




tay [8] 


* 


(INT64) 3 


+ tay [9] ; 





[000279] Once again, very few math operations are required compared to the matrix 
multiplication method. Thus, it becomes superfluous to use large amounts of memory to hold 
the transformation matrices. Also, since computation of the matrix values isn't performed at 
runtime, the transformations are extremely fast, 0(1). There is one disadvantage, hard-coding 
the transformations an create extremely large executable files; approximately 60 MB for both 
transformations through order 40. 

[000280] Plotting a wavefront using the Taylor monomials becomes trivial. First the user 
would define the Zernike polynomial in a vector format. Next the amplitudes would be 
transformed to Taylor amplitudes. Finally we can use the Taylor representation to plot the 
wavefront. 

c = (0,0, 0,0, 1,0) (158) 
a = fz*c 

n 

v{ x >y)= Yj a i xPix y Pi2 (159) 

where p is the vector holding the exponents for x and y. By using a class object we can write 

C++ code that looks almost identical to the mathematics above. Since all transforms will support 
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the same interface, many different transforms can be defined by using inheritance. If a class 
object is defined for transforms, we can then define what it means to multiply a transformation 
matrix by a vector of amplitudes by using operator overloading. The following C++ code does 
the same operation as the mathematics above 
[000281] 

CTransform *fm = new CTFZ (); // Define transform 

Wavefront W; // Define wavefront 

Poly c = { 0, 0, 0, 0, 1, 0 } // Define Z42 

Poly a = fz*c // Transform to Taylor 

for( int i=l; i<=n; i++ ) // Sum 

W[x] [y] = a[i] *pow(x, p[i][l]) *pow(y, p[i][2]) ; // Taylor components 

Delete fz; // Free transform 

[000282] First the code creates the class that transforms to the Taylor basis from the Zernike 
basis. Next we define a wavefront. The polynomial c is created to represent the eigenstate Z^2 
and then the Taylor amplitude vector a is created by multiplying c by the transform. Finally the 
wavefront is calculated in the Taylor representation. 

[000283] Using the said methods for code generation can yield hundreds of thousands of lines 
of code, depending on the order desired. Debugging that many lines of code is quite difficult, 
therefore, the use of mathematical validation is unavoidable. 

[000284] Several tests were done to show that the Zernike amplitudes were computed correctly 
without a single missed keystroke. For these tests, an independent programmer manually wrote 
C++ functions for the first five orders of the Zernike polynomials. So the production code and 
validation code were created by two different programmers. The first validation test was a point 



80 



Docket No. WFS.017 
--PATENT- 

by point test. This test compares the functional definitions of the Zernike polynomials through 
fifth order with the results of from using the Zernike array structure. The point by point test 
samples 10 M points within the unit circle and compared the values returned for both methods. 
There are no differences between the functional or the array representations. The next test uses 
the property of orthogonality to determine if the Zernike amplitudes were created correctly. The 
orthogonality of the Zernike polynomials, stated in equation 155, shows that if the product of 
two Zernike polynomials is integrated over the unit disc, the answer will be zero unless both 
polynomials are the same. The orthogonality test was a comprised of three parts. The first part 
tested the orthogonality of the explicit Zernike functions which were manually entered. The 
orthogonality of the explicit Zernike functions and the Zernike array structure was tested next. 
Finally the orthogonality of the Zernike array structure to itself was tested. In all three the results 
were difficult to interpret. The integration is preformed using the qquas routine. Because this 
is a numerical integration we get more and more error as the order increases. Because the first 
part of the test gave similar errors as the second and third, the non-zero terms off of the diagonal 
may be a result of the integration. Better methods of integration will have to be used to validate 
these conclusions. 

[000285] Proving that the transformation algorithms are correct can also be done in an 

automated manner. To determine accuracy of the algorithms we can create simple vectors for 

each Zernike eigenstate as well as the corresponding Taylor amplitude vector. A set of vectors 

through second order is shown in table 2. We transform the Zernike eigenstate into the Taylor 

basis and compare to the predicted value. To test the reverse transformation, the Taylor 

amplitudes are transformed to the Zernike basis and computed to the expected eigenvector. 
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Since this is all based on integer arithmetic, the differences are exactly zero. If there is a non- 
zero difference then the transformation algorithm has been programmed incorrectly. Even more 
interestingly, by using this method we can determine exactly where the problem occurred 
depending on which transformation yielded the non-zero result. 
[000286] Table 3: Simple Zernike Polynomials and the Equivalent Taylor Monomials 



n m 


Zernike 
Coefficients 


Taylor 
Coefficients 


0 0 


(1) 


(1) 


1 0 


(0,1,0) 


(0,1,0) 


1 1 


(0,0,1) 


(0,0,1) 


2 0 


(0,0,0,1,0,0) 


(0,0,0,02,0) 


2 1 


(0,0,0,0,1,0) 


(-1,0,0,2,0,2) 


2 2 


(0,0,0,0,0,1) 


(0,0,0,-1,0,1) 



[000287] Testing the precision of these affine transformations can be done by performing two 

simple tests repeatedly. These tests can give a good measure of how accurate these conversions 

are. To do these tests, create an array of type double to be transformed. The size of the array 

is (d + l)(d + 2)/2 where d is the order. Once the array has been created with the correct size, fill 

each element with a random number between -1 and 1. 

[000288] for(int ii = 0; ii < nTerms; ii++) 

initarray [ii] = ( (rand ()/ (double) RANDJVIAX) - 0.5) * 2; 
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[000289] With this array, assume that it contains Zemike amplitudes and perform the 
transformation to its Taylor equivalent. The intermediate array of Taylor amplitudes 
should then be converted back into its Zernike representation. 
[000290] intermediate = fz * initarray; 

finalarray = tz * intermediate; 
[000291] In theory, the finalarray should now contain the same values as our 
initarray. So, by doing this repeatedly for each other and with newly created random 
sequences each time, we can compare the sum of the squares of the differences between 
initarray and finalarray to determine the accuracy of our transformations. The second 
test is to repeat the above test doing the transformations in reverse. In other words, initarray 
should be initially assumed to be holding the Taylor coefficients. 

[000292] intermediate = tz * initarray; 

finalarray = fz * intermediate; 
[000293] The use of circle polynomials of Zernike in optical analysis provides valuable insight. 
Despite their role in optical analysis, computation using the Zernike polynomials is quite 
inefficient. We outlined how rational polynomials like those of Zernike, Legendre, Chebyshev 
and Laguerre can be handled as affine combinations of a Taylor monomial set. We demonstrated 
how calculations can be performed much more rapidly in the Taylor basis and how to use integer 
transformations to recover the exact amplitudes in the desired basis. We also discussed some 
basics of a C++ implementation. Finally testing and debugging a large volume of code requires 
mathematical validation. To insure the correctness of the code, validation was used to ensure the 

Zernike amplitudes were generated correctly and the transformations to and from Taylor worked 
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correctly. Even though Zernike polynomials can tell us a lot about an optical system they are not 
good for doing computation, therefore computation is done in Taylor space. 
[000294] While preferred embodiments are disclosed herein, many variations are possible 
which remain within the concept and scope of the invention. Such variations would become 
clearly to one of ordinary skill in the art after inspection of the specification, drawings and 
claims herein. The invention therefore is not to be restricted except within the spirit and scope of 
the appended claims. 
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